1 // Copyright (c) 2017, Oracle and/or its affiliates. All rights reserved.
2 //
3 // This program is free software; you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License, version 2.0,
5 // as published by the Free Software Foundation.
6 //
7 // This program is also distributed with certain software (including
8 // but not limited to OpenSSL) that is licensed under separate terms,
9 // as designated in a particular file or component or in included license
10 // documentation. The authors of MySQL hereby grant you an additional
11 // permission to link the program and your derivative works with the
12 // separately licensed software that they have included with MySQL.
13 //
14 // This program is distributed in the hope that it will be useful,
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 // GNU General Public License, version 2.0, for more details.
18 //
19 // You should have received a copy of the GNU General Public License
20 // along with this program; if not, write to the Free Software
21 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
22
23 /// @file
24 ///
25 /// This file implements the within functor and function.
26
27 #include <cmath> // std::isfinite
28 #include <limits>
29 #include <memory> // std::unique_ptr
30
31 #include <boost/geometry.hpp>
32
33 #include "sql/dd/types/spatial_reference_system.h" // dd::Spatial_reference_system
34 #include "sql/gis/box.h"
35 #include "sql/gis/box_traits.h"
36 #include "sql/gis/difference_functor.h"
37 #include "sql/gis/equals_functor.h"
38 #include "sql/gis/gc_utils.h"
39 #include "sql/gis/geometries.h"
40 #include "sql/gis/geometries_traits.h"
41 #include "sql/gis/intersects_functor.h"
42 #include "sql/gis/mbr_utils.h"
43 #include "sql/gis/relops.h"
44 #include "sql/gis/within_functor.h"
45 #include "sql/sql_exception_handler.h" // handle_gis_exception
46 #include "template_utils.h" // down_cast
47
48 namespace bg = boost::geometry;
49
50 namespace gis {
51
Within(double semi_major,double semi_minor)52 Within::Within(double semi_major, double semi_minor)
53 : m_semi_major(semi_major),
54 m_semi_minor(semi_minor),
55 m_geographic_pl_pa_strategy(
56 bg::srs::spheroid<double>(semi_major, semi_minor)),
57 m_geographic_ll_la_aa_strategy(
58 bg::srs::spheroid<double>(semi_major, semi_minor)) {}
59
operator ()(const Geometry * g1,const Geometry * g2) const60 bool Within::operator()(const Geometry *g1, const Geometry *g2) const {
61 return apply(*this, g1, g2);
62 }
63
operator ()(const Box * b1,const Box * b2) const64 bool Within::operator()(const Box *b1, const Box *b2) const {
65 DBUG_ASSERT(b1->coordinate_system() == b2->coordinate_system());
66 switch (b1->coordinate_system()) {
67 case Coordinate_system::kCartesian:
68 return eval(down_cast<const Cartesian_box *>(b1),
69 down_cast<const Cartesian_box *>(b2));
70 case Coordinate_system::kGeographic:
71 return eval(down_cast<const Geographic_box *>(b1),
72 down_cast<const Geographic_box *>(b2));
73 }
74
75 DBUG_ASSERT(false);
76 return false;
77 }
78
eval(const Geometry * g1,const Geometry * g2) const79 bool Within::eval(const Geometry *g1, const Geometry *g2) const {
80 // All parameter type combinations have been implemented.
81 DBUG_ASSERT(false);
82 throw not_implemented_exception::for_non_projected(*g1, *g2);
83 }
84
85 //////////////////////////////////////////////////////////////////////////////
86
87 // within(Cartesian_point, *)
88
eval(const Cartesian_point * g1,const Cartesian_point * g2) const89 bool Within::eval(const Cartesian_point *g1, const Cartesian_point *g2) const {
90 return bg::within(*g1, *g2);
91 }
92
eval(const Cartesian_point * g1,const Cartesian_linestring * g2) const93 bool Within::eval(const Cartesian_point *g1,
94 const Cartesian_linestring *g2) const {
95 return bg::within(*g1, *g2);
96 }
97
eval(const Cartesian_point * g1,const Cartesian_polygon * g2) const98 bool Within::eval(const Cartesian_point *g1,
99 const Cartesian_polygon *g2) const {
100 return bg::within(*g1, *g2);
101 }
102
eval(const Cartesian_point * g1,const Cartesian_geometrycollection * g2) const103 bool Within::eval(const Cartesian_point *g1,
104 const Cartesian_geometrycollection *g2) const {
105 // g1 must be within one of the elements of g2.
106 for (auto g : *g2)
107 if ((*this)(g1, g)) return true;
108 return false;
109 }
110
eval(const Cartesian_point * g1,const Cartesian_multipoint * g2) const111 bool Within::eval(const Cartesian_point *g1,
112 const Cartesian_multipoint *g2) const {
113 return bg::within(*g1, *g2);
114 }
115
eval(const Cartesian_point * g1,const Cartesian_multilinestring * g2) const116 bool Within::eval(const Cartesian_point *g1,
117 const Cartesian_multilinestring *g2) const {
118 return bg::within(*g1, *g2);
119 }
120
eval(const Cartesian_point * g1,const Cartesian_multipolygon * g2) const121 bool Within::eval(const Cartesian_point *g1,
122 const Cartesian_multipolygon *g2) const {
123 return bg::within(*g1, *g2);
124 }
125
126 //////////////////////////////////////////////////////////////////////////////
127
128 // within(Cartesian_linestring, *)
129
eval(const Cartesian_linestring *,const Cartesian_point *) const130 bool Within::eval(const Cartesian_linestring *, const Cartesian_point *) const {
131 // A linestring can never be within a point.
132 return false;
133 }
134
eval(const Cartesian_linestring * g1,const Cartesian_linestring * g2) const135 bool Within::eval(const Cartesian_linestring *g1,
136 const Cartesian_linestring *g2) const {
137 return bg::within(*g1, *g2);
138 }
139
eval(const Cartesian_linestring * g1,const Cartesian_polygon * g2) const140 bool Within::eval(const Cartesian_linestring *g1,
141 const Cartesian_polygon *g2) const {
142 return bg::within(*g1, *g2);
143 }
144
eval(const Cartesian_linestring * g1,const Cartesian_geometrycollection * g2) const145 bool Within::eval(const Cartesian_linestring *g1,
146 const Cartesian_geometrycollection *g2) const {
147 // For g1 to be within g2, no point of g1 may be in the exterior of g2 and at
148 // least one point of the interior of g1 has to be within the interior of g2.
149
150 std::unique_ptr<Multipoint> g2_mpt;
151 std::unique_ptr<Multilinestring> g2_mls;
152 std::unique_ptr<Multipolygon> g2_mpy;
153 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
154 &g2_mpy);
155 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
156
157 Difference difference(m_semi_major, m_semi_minor);
158 std::unique_ptr<Multilinestring> g1_diff_g2(new Cartesian_multilinestring());
159 g1_diff_g2.reset(down_cast<Cartesian_multilinestring *>(
160 difference(g1, down_cast<Cartesian_multilinestring *>(g2_mls.get()))));
161 g1_diff_g2.reset(down_cast<Cartesian_multilinestring *>(
162 difference(down_cast<Cartesian_multilinestring *>(g1_diff_g2.get()),
163 down_cast<Cartesian_multipolygon *>(g2_mpy.get()))));
164
165 boost::geometry::de9im::mask mask("T********");
166 return g1_diff_g2->empty() &&
167 (bg::relate(*g1, *down_cast<Cartesian_multilinestring *>(g2_mls.get()),
168 mask) ||
169 bg::relate(*g1, *down_cast<Cartesian_multipolygon *>(g2_mpy.get()),
170 mask));
171 }
172
eval(const Cartesian_linestring *,const Cartesian_multipoint *) const173 bool Within::eval(const Cartesian_linestring *,
174 const Cartesian_multipoint *) const {
175 // A linestring can never be within a multipoint.
176 return false;
177 }
178
eval(const Cartesian_linestring * g1,const Cartesian_multilinestring * g2) const179 bool Within::eval(const Cartesian_linestring *g1,
180 const Cartesian_multilinestring *g2) const {
181 return bg::within(*g1, *g2);
182 }
183
eval(const Cartesian_linestring * g1,const Cartesian_multipolygon * g2) const184 bool Within::eval(const Cartesian_linestring *g1,
185 const Cartesian_multipolygon *g2) const {
186 return bg::within(*g1, *g2);
187 }
188
189 //////////////////////////////////////////////////////////////////////////////
190
191 // within(Cartesian_polygon, *)
192
eval(const Cartesian_polygon *,const Cartesian_point *) const193 bool Within::eval(const Cartesian_polygon *, const Cartesian_point *) const {
194 // A polygon can never be within a point.
195 return false;
196 }
197
eval(const Cartesian_polygon *,const Cartesian_linestring *) const198 bool Within::eval(const Cartesian_polygon *,
199 const Cartesian_linestring *) const {
200 // A polygon can never be within a linestring.
201 return false;
202 }
203
eval(const Cartesian_polygon * g1,const Cartesian_polygon * g2) const204 bool Within::eval(const Cartesian_polygon *g1,
205 const Cartesian_polygon *g2) const {
206 return bg::within(*g1, *g2);
207 }
208
eval(const Cartesian_polygon * g1,const Cartesian_geometrycollection * g2) const209 bool Within::eval(const Cartesian_polygon *g1,
210 const Cartesian_geometrycollection *g2) const {
211 // Polygons can only be within other polygons or multipolygons, so we can
212 // ignore points and linestrings. g1 is within g2 if g1 is within the union
213 // multipolygon of g2.
214 std::unique_ptr<Multipoint> g2_mpt;
215 std::unique_ptr<Multilinestring> g2_mls;
216 std::unique_ptr<Multipolygon> g2_mpy;
217 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
218 &g2_mpy);
219 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
220 return eval(g1, down_cast<Cartesian_multipolygon *>(g2_mpy.get()));
221 }
222
eval(const Cartesian_polygon *,const Cartesian_multipoint *) const223 bool Within::eval(const Cartesian_polygon *,
224 const Cartesian_multipoint *) const {
225 // A polygon can never be within a multipoint.
226 return false;
227 }
228
eval(const Cartesian_polygon *,const Cartesian_multilinestring *) const229 bool Within::eval(const Cartesian_polygon *,
230 const Cartesian_multilinestring *) const {
231 // A polygon can never be within a multilinestring.
232 return false;
233 }
234
eval(const Cartesian_polygon * g1,const Cartesian_multipolygon * g2) const235 bool Within::eval(const Cartesian_polygon *g1,
236 const Cartesian_multipolygon *g2) const {
237 return bg::within(*g1, *g2);
238 }
239
240 //////////////////////////////////////////////////////////////////////////////
241
242 // within(Cartesian_geometrycollection, *)
243
eval(const Cartesian_geometrycollection * g1,const Cartesian_point * g2) const244 bool Within::eval(const Cartesian_geometrycollection *g1,
245 const Cartesian_point *g2) const {
246 Equals equals(m_semi_major, m_semi_minor);
247 return equals(g1, g2);
248 }
249
eval(const Cartesian_geometrycollection * g1,const Cartesian_linestring * g2) const250 bool Within::eval(const Cartesian_geometrycollection *g1,
251 const Cartesian_linestring *g2) const {
252 // g1 is within g2 if g1 contains only points and linestrings. One of the
253 // elements of g1 must be within g2, the rest must be covered by g2.
254 std::unique_ptr<Multipoint> g1_mpt;
255 std::unique_ptr<Multilinestring> g1_mls;
256 std::unique_ptr<Multipolygon> g1_mpy;
257 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
258 &g1_mpy);
259 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
260
261 if (!g1_mpy->empty()) return false;
262
263 if (eval(down_cast<Cartesian_multipoint *>(g1_mpt.get()), g2))
264 return g1_mls->empty() ||
265 bg::covered_by(*down_cast<Cartesian_multilinestring *>(g1_mls.get()),
266 *g2);
267 if (eval(down_cast<Cartesian_multilinestring *>(g1_mls.get()), g2)) {
268 for (auto &pt : *down_cast<Cartesian_multipoint *>(g1_mpt.get()))
269 if (!bg::covered_by(pt, *g2)) return false;
270 return true;
271 }
272 return false;
273 }
274
eval(const Cartesian_geometrycollection * g1,const Cartesian_polygon * g2) const275 bool Within::eval(const Cartesian_geometrycollection *g1,
276 const Cartesian_polygon *g2) const {
277 // At least one element of g1 has to be within g2. The rest have to be covered
278 // by g2.
279 std::unique_ptr<Multipoint> g1_mpt;
280 std::unique_ptr<Multilinestring> g1_mls;
281 std::unique_ptr<Multipolygon> g1_mpy;
282 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
283 &g1_mpy);
284 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
285
286 if (eval(down_cast<Cartesian_multipoint *>(g1_mpt.get()), g2)) {
287 return (g1_mls->empty() ||
288 bg::covered_by(
289 *down_cast<Cartesian_multilinestring *>(g1_mls.get()), *g2)) &&
290 (g1_mpy->empty() ||
291 bg::covered_by(*down_cast<Cartesian_multipolygon *>(g1_mpy.get()),
292 *g2));
293 }
294 if (eval(down_cast<Cartesian_multilinestring *>(g1_mls.get()), g2)) {
295 for (auto &pt : *down_cast<Cartesian_multipoint *>(g1_mpt.get()))
296 if (!bg::covered_by(pt, *g2)) return false;
297 return g1_mpy->empty() ||
298 bg::covered_by(*down_cast<Cartesian_multipolygon *>(g1_mpy.get()),
299 *g2);
300 }
301 if (eval(down_cast<Cartesian_multipolygon *>(g1_mpy.get()), g2)) {
302 for (auto &pt : *down_cast<Cartesian_multipoint *>(g1_mpt.get()))
303 if (!bg::covered_by(pt, *g2)) return false;
304 return g1_mls->empty() ||
305 bg::covered_by(*down_cast<Cartesian_multilinestring *>(g1_mls.get()),
306 *g2);
307 }
308 return false;
309 }
310
eval(const Cartesian_geometrycollection * g1,const Cartesian_geometrycollection * g2) const311 bool Within::eval(const Cartesian_geometrycollection *g1,
312 const Cartesian_geometrycollection *g2) const {
313 // At least one element of g1 has to be within g2. The rest have to be covered
314 // by g2.
315 std::unique_ptr<Multipoint> g1_mpt;
316 std::unique_ptr<Multilinestring> g1_mls;
317 std::unique_ptr<Multipolygon> g1_mpy;
318 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
319 &g1_mpy);
320 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
321
322 std::unique_ptr<Multipoint> g2_mpt;
323 std::unique_ptr<Multilinestring> g2_mls;
324 std::unique_ptr<Multipolygon> g2_mpy;
325 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
326 &g2_mpy);
327 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
328
329 // Check that no part of g1 is in the exterior of g2.
330 Difference difference(m_semi_major, m_semi_minor);
331 std::unique_ptr<Cartesian_multipoint> g1_mpt_diff_g2(
332 new Cartesian_multipoint());
333 g1_mpt_diff_g2.reset(down_cast<Cartesian_multipoint *>(
334 difference(down_cast<Cartesian_multipoint *>(g1_mpt.get()),
335 down_cast<Cartesian_multipoint *>(g2_mpt.get()))));
336 g1_mpt_diff_g2.reset(down_cast<Cartesian_multipoint *>(
337 difference(g1_mpt_diff_g2.get(),
338 down_cast<Cartesian_multilinestring *>(g2_mls.get()))));
339 g1_mpt_diff_g2.reset(down_cast<Cartesian_multipoint *>(
340 difference(g1_mpt_diff_g2.get(),
341 down_cast<Cartesian_multipolygon *>(g2_mpy.get()))));
342 if (!g1_mpt_diff_g2->empty()) return false;
343
344 std::unique_ptr<Cartesian_multilinestring> g1_mls_diff_g2(
345 new Cartesian_multilinestring());
346 g1_mls_diff_g2.reset(down_cast<Cartesian_multilinestring *>(
347 difference(down_cast<Cartesian_multilinestring *>(g1_mls.get()),
348 down_cast<Cartesian_multilinestring *>(g2_mls.get()))));
349 g1_mls_diff_g2.reset(down_cast<Cartesian_multilinestring *>(
350 difference(g1_mls_diff_g2.get(),
351 down_cast<Cartesian_multipolygon *>(g2_mpy.get()))));
352 if (!g1_mls_diff_g2->empty()) return false;
353
354 std::unique_ptr<Cartesian_multipolygon> g1_mpy_diff_g2(
355 new Cartesian_multipolygon());
356 g1_mpy_diff_g2.reset(down_cast<Cartesian_multipolygon *>(
357 difference(down_cast<Cartesian_multipolygon *>(g1_mpy.get()),
358 down_cast<Cartesian_multipolygon *>(g2_mpy.get()))));
359 if (!g1_mpy_diff_g2->empty()) return false;
360
361 // Check that the interiors of g1 and g2 have at least one point in common.
362 boost::geometry::de9im::mask mask("T********");
363 return eval(down_cast<Cartesian_multipoint *>(g1_mpt.get()), g2) ||
364 bg::relate(*down_cast<Cartesian_multilinestring *>(g1_mls.get()),
365 *down_cast<Cartesian_multilinestring *>(g2_mls.get()),
366 mask) ||
367 bg::relate(*down_cast<Cartesian_multilinestring *>(g1_mls.get()),
368 *down_cast<Cartesian_multipolygon *>(g2_mpy.get()), mask) ||
369 bg::relate(*down_cast<Cartesian_multipolygon *>(g1_mpy.get()),
370 *down_cast<Cartesian_multipolygon *>(g2_mpy.get()), mask);
371 }
372
eval(const Cartesian_geometrycollection * g1,const Cartesian_multipoint * g2) const373 bool Within::eval(const Cartesian_geometrycollection *g1,
374 const Cartesian_multipoint *g2) const {
375 // g1 is within g2 if g1 contains only points and those points are within g2.
376 std::unique_ptr<Multipoint> g1_mpt;
377 std::unique_ptr<Multilinestring> g1_mls;
378 std::unique_ptr<Multipolygon> g1_mpy;
379 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
380 &g1_mpy);
381 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
382 return g1_mls->empty() && g1_mpy->empty() &&
383 eval(down_cast<Cartesian_multipoint *>(g1_mpt.get()), g2);
384 }
385
eval(const Cartesian_geometrycollection * g1,const Cartesian_multilinestring * g2) const386 bool Within::eval(const Cartesian_geometrycollection *g1,
387 const Cartesian_multilinestring *g2) const {
388 // g1 is within g2 if g1 contains only points and linestrings. One of the
389 // elements of g1 must be within g2, the rest must be covered by g2.
390 std::unique_ptr<Multipoint> g1_mpt;
391 std::unique_ptr<Multilinestring> g1_mls;
392 std::unique_ptr<Multipolygon> g1_mpy;
393 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
394 &g1_mpy);
395 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
396
397 if (!g1_mpy->empty()) return false;
398
399 if (eval(down_cast<Cartesian_multipoint *>(g1_mpt.get()), g2))
400 return g1_mls->empty() ||
401 bg::covered_by(*down_cast<Cartesian_multilinestring *>(g1_mls.get()),
402 *g2);
403 if (eval(down_cast<Cartesian_multilinestring *>(g1_mls.get()), g2)) {
404 for (auto &pt : *down_cast<Cartesian_multipoint *>(g1_mpt.get()))
405 if (!bg::covered_by(pt, *g2)) return false;
406 return true;
407 }
408 return false;
409 }
410
eval(const Cartesian_geometrycollection * g1,const Cartesian_multipolygon * g2) const411 bool Within::eval(const Cartesian_geometrycollection *g1,
412 const Cartesian_multipolygon *g2) const {
413 // At least one element of g1 has to be within g2. The rest have to be covered
414 // by g2.
415 std::unique_ptr<Multipoint> g1_mpt;
416 std::unique_ptr<Multilinestring> g1_mls;
417 std::unique_ptr<Multipolygon> g1_mpy;
418 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
419 &g1_mpy);
420 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
421
422 if (eval(down_cast<Cartesian_multipoint *>(g1_mpt.get()), g2)) {
423 return (g1_mls->empty() ||
424 bg::covered_by(
425 *down_cast<Cartesian_multilinestring *>(g1_mls.get()), *g2)) &&
426 (g1_mpy->empty() ||
427 bg::covered_by(*down_cast<Cartesian_multipolygon *>(g1_mpy.get()),
428 *g2));
429 }
430 if (eval(down_cast<Cartesian_multilinestring *>(g1_mls.get()), g2)) {
431 for (auto &pt : *down_cast<Cartesian_multipoint *>(g1_mpt.get()))
432 if (!bg::covered_by(pt, *g2)) return false;
433 return g1_mpy->empty() ||
434 bg::covered_by(*down_cast<Cartesian_multipolygon *>(g1_mpy.get()),
435 *g2);
436 }
437 if (eval(down_cast<Cartesian_multipolygon *>(g1_mpy.get()), g2)) {
438 for (auto &pt : *down_cast<Cartesian_multipoint *>(g1_mpt.get()))
439 if (!bg::covered_by(pt, *g2)) return false;
440 return g1_mls->empty() ||
441 bg::covered_by(*down_cast<Cartesian_multilinestring *>(g1_mls.get()),
442 *g2);
443 }
444 return false;
445 }
446
447 //////////////////////////////////////////////////////////////////////////////
448
449 // within(Cartesian_multipoint, *)
450
eval(const Cartesian_multipoint * g1,const Cartesian_point * g2) const451 bool Within::eval(const Cartesian_multipoint *g1,
452 const Cartesian_point *g2) const {
453 Equals equals(m_semi_major, m_semi_minor);
454 return equals(g1, g2);
455 }
456
eval(const Cartesian_multipoint * g1,const Cartesian_linestring * g2) const457 bool Within::eval(const Cartesian_multipoint *g1,
458 const Cartesian_linestring *g2) const {
459 // At least one point in g1 must be within g2. The rest has to intersect g2.
460 bool within = false;
461 bool intersects = false;
462 for (auto &pt : *g1) {
463 if (!within) {
464 within = bg::within(pt, *g2);
465 if (!within)
466 intersects = bg::intersects(pt, *g2);
467 else
468 intersects = true;
469 } else {
470 intersects = bg::intersects(pt, *g2);
471 }
472 if (!intersects) break;
473 }
474 return (within && intersects);
475 }
476
eval(const Cartesian_multipoint * g1,const Cartesian_polygon * g2) const477 bool Within::eval(const Cartesian_multipoint *g1,
478 const Cartesian_polygon *g2) const {
479 // At least one point in g1 must be within g2. The rest has to intersect g2.
480 bool within = false;
481 bool intersects = false;
482 for (auto &pt : *g1) {
483 if (!within) {
484 within = bg::within(pt, *g2);
485 if (!within)
486 intersects = bg::intersects(pt, *g2);
487 else
488 intersects = true;
489 } else {
490 intersects = bg::intersects(pt, *g2);
491 }
492 if (!intersects) break;
493 }
494 return (within && intersects);
495 }
496
eval(const Cartesian_multipoint * g1,const Cartesian_geometrycollection * g2) const497 bool Within::eval(const Cartesian_multipoint *g1,
498 const Cartesian_geometrycollection *g2) const {
499 // At least one point in g1 must be within g2. The rest has to intersect g2.
500 Intersects intersects_func(m_semi_major, m_semi_minor);
501 bool within = false;
502 bool intersects = false;
503 for (auto &pt : *g1) {
504 if (!within) {
505 within = eval(&pt, g2);
506 if (!within)
507 intersects = intersects_func(&pt, g2);
508 else
509 intersects = true;
510 } else {
511 intersects = intersects_func(&pt, g2);
512 }
513 if (!intersects) break;
514 }
515 return (within && intersects);
516 }
517
eval(const Cartesian_multipoint * g1,const Cartesian_multipoint * g2) const518 bool Within::eval(const Cartesian_multipoint *g1,
519 const Cartesian_multipoint *g2) const {
520 for (auto &pt : *g1)
521 if (!bg::within(pt, *g2)) return false;
522 return true;
523 }
524
eval(const Cartesian_multipoint * g1,const Cartesian_multilinestring * g2) const525 bool Within::eval(const Cartesian_multipoint *g1,
526 const Cartesian_multilinestring *g2) const {
527 // At least one point in g1 must be within g2. The rest has to intersect g2.
528 bool within = false;
529 bool intersects = false;
530 for (auto &pt : *g1) {
531 if (!within) {
532 within = bg::within(pt, *g2);
533 if (!within)
534 intersects = bg::intersects(pt, *g2);
535 else
536 intersects = true;
537 } else {
538 intersects = bg::intersects(pt, *g2);
539 }
540 if (!intersects) break;
541 }
542 return (within && intersects);
543 }
544
eval(const Cartesian_multipoint * g1,const Cartesian_multipolygon * g2) const545 bool Within::eval(const Cartesian_multipoint *g1,
546 const Cartesian_multipolygon *g2) const {
547 // At least one point in g1 must be within g2. The rest has to intersect g2.
548 bool within = false;
549 bool intersects = false;
550 for (auto &pt : *g1) {
551 if (!within) {
552 within = bg::within(pt, *g2);
553 if (!within)
554 intersects = bg::intersects(pt, *g2);
555 else
556 intersects = true;
557 } else {
558 intersects = bg::intersects(pt, *g2);
559 }
560 if (!intersects) break;
561 }
562 return (within && intersects);
563 }
564
565 //////////////////////////////////////////////////////////////////////////////
566
567 // within(Cartesian_multilinestring, *)
568
eval(const Cartesian_multilinestring *,const Cartesian_point *) const569 bool Within::eval(const Cartesian_multilinestring *,
570 const Cartesian_point *) const {
571 // A multilinestring can never be within a point.
572 return false;
573 }
574
eval(const Cartesian_multilinestring * g1,const Cartesian_linestring * g2) const575 bool Within::eval(const Cartesian_multilinestring *g1,
576 const Cartesian_linestring *g2) const {
577 return bg::within(*g1, *g2);
578 }
579
eval(const Cartesian_multilinestring * g1,const Cartesian_polygon * g2) const580 bool Within::eval(const Cartesian_multilinestring *g1,
581 const Cartesian_polygon *g2) const {
582 return bg::within(*g1, *g2);
583 }
584
eval(const Cartesian_multilinestring * g1,const Cartesian_geometrycollection * g2) const585 bool Within::eval(const Cartesian_multilinestring *g1,
586 const Cartesian_geometrycollection *g2) const {
587 // For g1 to be within g2, no point of g1 may be in the exterior of g2 and at
588 // least one point of the interior of g1 has to be within the interior of g2.
589
590 std::unique_ptr<Multipoint> g2_mpt;
591 std::unique_ptr<Multilinestring> g2_mls;
592 std::unique_ptr<Multipolygon> g2_mpy;
593 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
594 &g2_mpy);
595 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
596
597 Difference difference(m_semi_major, m_semi_minor);
598 std::unique_ptr<Cartesian_multilinestring> g1_diff_g2(
599 new Cartesian_multilinestring());
600 g1_diff_g2.reset(down_cast<Cartesian_multilinestring *>(
601 difference(g1, down_cast<Cartesian_multilinestring *>(g2_mls.get()))));
602 g1_diff_g2.reset(down_cast<Cartesian_multilinestring *>(difference(
603 g1_diff_g2.get(), down_cast<Cartesian_multipolygon *>(g2_mpy.get()))));
604
605 boost::geometry::de9im::mask mask("T********");
606 return g1_diff_g2->empty() &&
607 (bg::relate(*g1, *down_cast<Cartesian_multilinestring *>(g2_mls.get()),
608 mask) ||
609 bg::relate(*g1, *down_cast<Cartesian_multipolygon *>(g2_mpy.get()),
610 mask));
611 }
612
eval(const Cartesian_multilinestring *,const Cartesian_multipoint *) const613 bool Within::eval(const Cartesian_multilinestring *,
614 const Cartesian_multipoint *) const {
615 // A multilinestring can never be within a multipoint.
616 return false;
617 }
618
eval(const Cartesian_multilinestring * g1,const Cartesian_multilinestring * g2) const619 bool Within::eval(const Cartesian_multilinestring *g1,
620 const Cartesian_multilinestring *g2) const {
621 return bg::within(*g1, *g2);
622 }
623
eval(const Cartesian_multilinestring * g1,const Cartesian_multipolygon * g2) const624 bool Within::eval(const Cartesian_multilinestring *g1,
625 const Cartesian_multipolygon *g2) const {
626 return bg::within(*g1, *g2);
627 }
628
629 //////////////////////////////////////////////////////////////////////////////
630
631 // within(Cartesian_multipolygon, *)
632
eval(const Cartesian_multipolygon *,const Cartesian_point *) const633 bool Within::eval(const Cartesian_multipolygon *,
634 const Cartesian_point *) const {
635 // A multipolygon can never be within a point.
636 return false;
637 }
638
eval(const Cartesian_multipolygon *,const Cartesian_linestring *) const639 bool Within::eval(const Cartesian_multipolygon *,
640 const Cartesian_linestring *) const {
641 // A multipolygon can never be within a linestring.
642 return false;
643 }
644
eval(const Cartesian_multipolygon * g1,const Cartesian_polygon * g2) const645 bool Within::eval(const Cartesian_multipolygon *g1,
646 const Cartesian_polygon *g2) const {
647 return bg::within(*g1, *g2);
648 }
649
eval(const Cartesian_multipolygon * g1,const Cartesian_geometrycollection * g2) const650 bool Within::eval(const Cartesian_multipolygon *g1,
651 const Cartesian_geometrycollection *g2) const {
652 // A multipolygon may not be within the points and linestrings of g2, so the
653 // only way a multipolygon is within a geometrycollectin, is if it's within
654 // the union multipolygon of the geometrycollection.
655 std::unique_ptr<Multipoint> g2_mpt;
656 std::unique_ptr<Multilinestring> g2_mls;
657 std::unique_ptr<Multipolygon> g2_mpy;
658 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
659 &g2_mpy);
660 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
661 return eval(g1, down_cast<Cartesian_multipolygon *>(g2_mpy.get()));
662 }
663
eval(const Cartesian_multipolygon *,const Cartesian_multipoint *) const664 bool Within::eval(const Cartesian_multipolygon *,
665 const Cartesian_multipoint *) const {
666 // A multipolygon can never be within a multipoint.
667 return false;
668 }
669
eval(const Cartesian_multipolygon *,const Cartesian_multilinestring *) const670 bool Within::eval(const Cartesian_multipolygon *,
671 const Cartesian_multilinestring *) const {
672 // A multipolygon can never be within a multilinestring.
673 return false;
674 }
675
eval(const Cartesian_multipolygon * g1,const Cartesian_multipolygon * g2) const676 bool Within::eval(const Cartesian_multipolygon *g1,
677 const Cartesian_multipolygon *g2) const {
678 return bg::within(*g1, *g2);
679 }
680
681 //////////////////////////////////////////////////////////////////////////////
682
683 // within(Geographic_point, *)
684
eval(const Geographic_point * g1,const Geographic_point * g2) const685 bool Within::eval(const Geographic_point *g1,
686 const Geographic_point *g2) const {
687 // Default strategy is OK. P/P computations do not depend on shape of
688 // ellipsoid.
689 return bg::within(*g1, *g2);
690 }
691
eval(const Geographic_point * g1,const Geographic_linestring * g2) const692 bool Within::eval(const Geographic_point *g1,
693 const Geographic_linestring *g2) const {
694 return bg::within(*g1, *g2, m_geographic_pl_pa_strategy);
695 }
696
eval(const Geographic_point * g1,const Geographic_polygon * g2) const697 bool Within::eval(const Geographic_point *g1,
698 const Geographic_polygon *g2) const {
699 return bg::within(*g1, *g2, m_geographic_pl_pa_strategy);
700 }
701
eval(const Geographic_point * g1,const Geographic_geometrycollection * g2) const702 bool Within::eval(const Geographic_point *g1,
703 const Geographic_geometrycollection *g2) const {
704 // g1 must be within one of the elements of g2.
705 for (auto g : *g2)
706 if ((*this)(g1, g)) return true;
707 return false;
708 }
709
eval(const Geographic_point * g1,const Geographic_multipoint * g2) const710 bool Within::eval(const Geographic_point *g1,
711 const Geographic_multipoint *g2) const {
712 // Default strategy is OK. P/P computations do not depend on shape of
713 // ellipsoid.
714 return bg::within(*g1, *g2);
715 }
716
eval(const Geographic_point * g1,const Geographic_multilinestring * g2) const717 bool Within::eval(const Geographic_point *g1,
718 const Geographic_multilinestring *g2) const {
719 return bg::within(*g1, *g2, m_geographic_pl_pa_strategy);
720 }
721
eval(const Geographic_point * g1,const Geographic_multipolygon * g2) const722 bool Within::eval(const Geographic_point *g1,
723 const Geographic_multipolygon *g2) const {
724 return bg::within(*g1, *g2, m_geographic_pl_pa_strategy);
725 }
726
727 //////////////////////////////////////////////////////////////////////////////
728
729 // within(Geographic_linestring, *)
730
eval(const Geographic_linestring *,const Geographic_point *) const731 bool Within::eval(const Geographic_linestring *,
732 const Geographic_point *) const {
733 // A linestring can never be within a point.
734 return false;
735 }
736
eval(const Geographic_linestring * g1,const Geographic_linestring * g2) const737 bool Within::eval(const Geographic_linestring *g1,
738 const Geographic_linestring *g2) const {
739 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
740 }
741
eval(const Geographic_linestring * g1,const Geographic_polygon * g2) const742 bool Within::eval(const Geographic_linestring *g1,
743 const Geographic_polygon *g2) const {
744 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
745 }
746
eval(const Geographic_linestring * g1,const Geographic_geometrycollection * g2) const747 bool Within::eval(const Geographic_linestring *g1,
748 const Geographic_geometrycollection *g2) const {
749 // For g1 to be within g2, no point of g1 may be in the exterior of g2 and at
750 // least one point of the interior of g1 has to be within the interior of g2.
751
752 std::unique_ptr<Multipoint> g2_mpt;
753 std::unique_ptr<Multilinestring> g2_mls;
754 std::unique_ptr<Multipolygon> g2_mpy;
755 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
756 &g2_mpy);
757 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
758
759 Difference difference(m_semi_major, m_semi_minor);
760 std::unique_ptr<Multilinestring> g1_diff_g2(new Geographic_multilinestring());
761 g1_diff_g2.reset(down_cast<Geographic_multilinestring *>(
762 difference(g1, down_cast<Geographic_multilinestring *>(g2_mls.get()))));
763 g1_diff_g2.reset(down_cast<Geographic_multilinestring *>(
764 difference(down_cast<Geographic_multilinestring *>(g1_diff_g2.get()),
765 down_cast<Geographic_multipolygon *>(g2_mpy.get()))));
766
767 boost::geometry::de9im::mask mask("T********");
768 return g1_diff_g2->empty() &&
769 (bg::relate(*g1,
770 *down_cast<Geographic_multilinestring *>(g2_mls.get()),
771 mask, m_geographic_ll_la_aa_strategy) ||
772 bg::relate(*g1, *down_cast<Geographic_multipolygon *>(g2_mpy.get()),
773 mask, m_geographic_ll_la_aa_strategy));
774 }
775
eval(const Geographic_linestring *,const Geographic_multipoint *) const776 bool Within::eval(const Geographic_linestring *,
777 const Geographic_multipoint *) const {
778 // A linestring can never be within a multipoint.
779 return false;
780 }
781
eval(const Geographic_linestring * g1,const Geographic_multilinestring * g2) const782 bool Within::eval(const Geographic_linestring *g1,
783 const Geographic_multilinestring *g2) const {
784 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
785 }
786
eval(const Geographic_linestring * g1,const Geographic_multipolygon * g2) const787 bool Within::eval(const Geographic_linestring *g1,
788 const Geographic_multipolygon *g2) const {
789 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
790 }
791
792 //////////////////////////////////////////////////////////////////////////////
793
794 // within(Geographic_polygon, *)
795
eval(const Geographic_polygon *,const Geographic_point *) const796 bool Within::eval(const Geographic_polygon *, const Geographic_point *) const {
797 // A polygon can never be within a point.
798 return false;
799 }
800
eval(const Geographic_polygon *,const Geographic_linestring *) const801 bool Within::eval(const Geographic_polygon *,
802 const Geographic_linestring *) const {
803 // A polygon can never be within a linestring.
804 return false;
805 }
806
eval(const Geographic_polygon * g1,const Geographic_polygon * g2) const807 bool Within::eval(const Geographic_polygon *g1,
808 const Geographic_polygon *g2) const {
809 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
810 }
811
eval(const Geographic_polygon * g1,const Geographic_geometrycollection * g2) const812 bool Within::eval(const Geographic_polygon *g1,
813 const Geographic_geometrycollection *g2) const {
814 // Polygons can only be within other polygons or multipolygons, so we can
815 // ignore points and linestrings. g1 is within g2 if g1 is within the union
816 // multipolygon of g2.
817 std::unique_ptr<Multipoint> g2_mpt;
818 std::unique_ptr<Multilinestring> g2_mls;
819 std::unique_ptr<Multipolygon> g2_mpy;
820 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
821 &g2_mpy);
822 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
823 return eval(g1, down_cast<Geographic_multipolygon *>(g2_mpy.get()));
824 }
825
eval(const Geographic_polygon *,const Geographic_multipoint *) const826 bool Within::eval(const Geographic_polygon *,
827 const Geographic_multipoint *) const {
828 // A polygon can never be within a multipoint.
829 return false;
830 }
831
eval(const Geographic_polygon *,const Geographic_multilinestring *) const832 bool Within::eval(const Geographic_polygon *,
833 const Geographic_multilinestring *) const {
834 // A polygon can never be within a multilinestring.
835 return false;
836 }
837
eval(const Geographic_polygon * g1,const Geographic_multipolygon * g2) const838 bool Within::eval(const Geographic_polygon *g1,
839 const Geographic_multipolygon *g2) const {
840 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
841 }
842
843 //////////////////////////////////////////////////////////////////////////////
844
845 // within(Geographic_geometrycollection, *)
846
eval(const Geographic_geometrycollection * g1,const Geographic_point * g2) const847 bool Within::eval(const Geographic_geometrycollection *g1,
848 const Geographic_point *g2) const {
849 Equals equals(m_semi_major, m_semi_minor);
850 return equals(g1, g2);
851 }
852
eval(const Geographic_geometrycollection * g1,const Geographic_linestring * g2) const853 bool Within::eval(const Geographic_geometrycollection *g1,
854 const Geographic_linestring *g2) const {
855 // g1 is within g2 if g1 contains only points and linestrings. One of the
856 // elements of g1 must be within g2, the rest must be covered by g2.
857 std::unique_ptr<Multipoint> g1_mpt;
858 std::unique_ptr<Multilinestring> g1_mls;
859 std::unique_ptr<Multipolygon> g1_mpy;
860 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
861 &g1_mpy);
862 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
863
864 if (!g1_mpy->empty()) return false;
865
866 if (eval(down_cast<Geographic_multipoint *>(g1_mpt.get()), g2))
867 return g1_mls->empty() ||
868 bg::covered_by(
869 *down_cast<Geographic_multilinestring *>(g1_mls.get()), *g2,
870 m_geographic_ll_la_aa_strategy);
871 if (eval(down_cast<Geographic_multilinestring *>(g1_mls.get()), g2)) {
872 for (auto &pt : *down_cast<Geographic_multipoint *>(g1_mpt.get()))
873 if (!bg::covered_by(pt, *g2, m_geographic_pl_pa_strategy)) return false;
874 return true;
875 }
876 return false;
877 }
878
eval(const Geographic_geometrycollection * g1,const Geographic_polygon * g2) const879 bool Within::eval(const Geographic_geometrycollection *g1,
880 const Geographic_polygon *g2) const {
881 // At least one element of g1 has to be within g2. The rest have to be covered
882 // by g2.
883 std::unique_ptr<Multipoint> g1_mpt;
884 std::unique_ptr<Multilinestring> g1_mls;
885 std::unique_ptr<Multipolygon> g1_mpy;
886 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
887 &g1_mpy);
888 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
889
890 if (eval(down_cast<Geographic_multipoint *>(g1_mpt.get()), g2)) {
891 return (g1_mls->empty() ||
892 bg::covered_by(
893 *down_cast<Geographic_multilinestring *>(g1_mls.get()), *g2,
894 m_geographic_ll_la_aa_strategy)) &&
895 (g1_mpy->empty() ||
896 bg::covered_by(*down_cast<Geographic_multipolygon *>(g1_mpy.get()),
897 *g2, m_geographic_ll_la_aa_strategy));
898 }
899 if (eval(down_cast<Geographic_multilinestring *>(g1_mls.get()), g2)) {
900 for (auto &pt : *down_cast<Geographic_multipoint *>(g1_mpt.get()))
901 if (!bg::covered_by(pt, *g2, m_geographic_pl_pa_strategy)) return false;
902 return g1_mpy->empty() ||
903 bg::covered_by(*down_cast<Geographic_multipolygon *>(g1_mpy.get()),
904 *g2, m_geographic_ll_la_aa_strategy);
905 }
906 if (eval(down_cast<Geographic_multipolygon *>(g1_mpy.get()), g2)) {
907 for (auto &pt : *down_cast<Geographic_multipoint *>(g1_mpt.get()))
908 if (!bg::covered_by(pt, *g2, m_geographic_pl_pa_strategy)) return false;
909 return g1_mls->empty() ||
910 bg::covered_by(
911 *down_cast<Geographic_multilinestring *>(g1_mls.get()), *g2,
912 m_geographic_ll_la_aa_strategy);
913 }
914 return false;
915 }
916
eval(const Geographic_geometrycollection * g1,const Geographic_geometrycollection * g2) const917 bool Within::eval(const Geographic_geometrycollection *g1,
918 const Geographic_geometrycollection *g2) const {
919 // At least one element of g1 has to be within g2. The rest have to be covered
920 // by g2.
921 std::unique_ptr<Multipoint> g1_mpt;
922 std::unique_ptr<Multilinestring> g1_mls;
923 std::unique_ptr<Multipolygon> g1_mpy;
924 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
925 &g1_mpy);
926 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
927
928 std::unique_ptr<Multipoint> g2_mpt;
929 std::unique_ptr<Multilinestring> g2_mls;
930 std::unique_ptr<Multipolygon> g2_mpy;
931 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
932 &g2_mpy);
933 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
934
935 // Check that no part of g1 is in the exterior of g2.
936 Difference difference(m_semi_major, m_semi_minor);
937 std::unique_ptr<Geographic_multipoint> g1_mpt_diff_g2(
938 new Geographic_multipoint());
939 g1_mpt_diff_g2.reset(down_cast<Geographic_multipoint *>(
940 difference(down_cast<Geographic_multipoint *>(g1_mpt.get()),
941 down_cast<Geographic_multipoint *>(g2_mpt.get()))));
942 g1_mpt_diff_g2.reset(down_cast<Geographic_multipoint *>(
943 difference(g1_mpt_diff_g2.get(),
944 down_cast<Geographic_multilinestring *>(g2_mls.get()))));
945 g1_mpt_diff_g2.reset(down_cast<Geographic_multipoint *>(
946 difference(g1_mpt_diff_g2.get(),
947 down_cast<Geographic_multipolygon *>(g2_mpy.get()))));
948 if (!g1_mpt_diff_g2->empty()) return false;
949
950 std::unique_ptr<Geographic_multilinestring> g1_mls_diff_g2(
951 new Geographic_multilinestring());
952 g1_mls_diff_g2.reset(down_cast<Geographic_multilinestring *>(
953 difference(down_cast<Geographic_multilinestring *>(g1_mls.get()),
954 down_cast<Geographic_multilinestring *>(g2_mls.get()))));
955 g1_mls_diff_g2.reset(down_cast<Geographic_multilinestring *>(
956 difference(g1_mls_diff_g2.get(),
957 down_cast<Geographic_multipolygon *>(g2_mpy.get()))));
958 if (!g1_mls_diff_g2->empty()) return false;
959
960 std::unique_ptr<Geographic_multipolygon> g1_mpy_diff_g2(
961 new Geographic_multipolygon());
962 g1_mpy_diff_g2.reset(down_cast<Geographic_multipolygon *>(
963 difference(down_cast<Geographic_multipolygon *>(g1_mpy.get()),
964 down_cast<Geographic_multipolygon *>(g2_mpy.get()))));
965 if (!g1_mpy_diff_g2->empty()) return false;
966
967 // Check that the interiors of g1 and g2 have at least one point in common.
968 boost::geometry::de9im::mask mask("T********");
969 return eval(down_cast<Geographic_multipoint *>(g1_mpt.get()), g2) ||
970 bg::relate(*down_cast<Geographic_multilinestring *>(g1_mls.get()),
971 *down_cast<Geographic_multilinestring *>(g2_mls.get()),
972 mask, m_geographic_ll_la_aa_strategy) ||
973 bg::relate(*down_cast<Geographic_multilinestring *>(g1_mls.get()),
974 *down_cast<Geographic_multipolygon *>(g2_mpy.get()), mask,
975 m_geographic_ll_la_aa_strategy) ||
976 bg::relate(*down_cast<Geographic_multipolygon *>(g1_mpy.get()),
977 *down_cast<Geographic_multipolygon *>(g2_mpy.get()), mask,
978 m_geographic_ll_la_aa_strategy);
979 }
980
eval(const Geographic_geometrycollection * g1,const Geographic_multipoint * g2) const981 bool Within::eval(const Geographic_geometrycollection *g1,
982 const Geographic_multipoint *g2) const {
983 // g1 is within g2 if g1 contains only points and those points are within g2.
984 std::unique_ptr<Multipoint> g1_mpt;
985 std::unique_ptr<Multilinestring> g1_mls;
986 std::unique_ptr<Multipolygon> g1_mpy;
987 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
988 &g1_mpy);
989 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
990 return g1_mls->empty() && g1_mpy->empty() &&
991 eval(down_cast<Geographic_multipoint *>(g1_mpt.get()), g2);
992 }
993
eval(const Geographic_geometrycollection * g1,const Geographic_multilinestring * g2) const994 bool Within::eval(const Geographic_geometrycollection *g1,
995 const Geographic_multilinestring *g2) const {
996 // g1 is within g2 if g1 contains only points and linestrings. One of the
997 // elements of g1 must be within g2, the rest must be covered by g2.
998 std::unique_ptr<Multipoint> g1_mpt;
999 std::unique_ptr<Multilinestring> g1_mls;
1000 std::unique_ptr<Multipolygon> g1_mpy;
1001 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
1002 &g1_mpy);
1003 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
1004
1005 if (!g1_mpy->empty()) return false;
1006
1007 if (eval(down_cast<Geographic_multipoint *>(g1_mpt.get()), g2))
1008 return g1_mls->empty() ||
1009 bg::covered_by(
1010 *down_cast<Geographic_multilinestring *>(g1_mls.get()), *g2,
1011 m_geographic_ll_la_aa_strategy);
1012 if (eval(down_cast<Geographic_multilinestring *>(g1_mls.get()), g2)) {
1013 for (auto &pt : *down_cast<Geographic_multipoint *>(g1_mpt.get()))
1014 if (!bg::covered_by(pt, *g2, m_geographic_pl_pa_strategy)) return false;
1015 return true;
1016 }
1017 return false;
1018 }
1019
eval(const Geographic_geometrycollection * g1,const Geographic_multipolygon * g2) const1020 bool Within::eval(const Geographic_geometrycollection *g1,
1021 const Geographic_multipolygon *g2) const {
1022 // At least one element of g1 has to be within g2. The rest have to be covered
1023 // by g2.
1024 std::unique_ptr<Multipoint> g1_mpt;
1025 std::unique_ptr<Multilinestring> g1_mls;
1026 std::unique_ptr<Multipolygon> g1_mpy;
1027 split_gc(down_cast<const Geometrycollection *>(g1), &g1_mpt, &g1_mls,
1028 &g1_mpy);
1029 gc_union(m_semi_major, m_semi_minor, &g1_mpt, &g1_mls, &g1_mpy);
1030
1031 if (eval(down_cast<Geographic_multipoint *>(g1_mpt.get()), g2)) {
1032 return (g1_mls->empty() ||
1033 bg::covered_by(
1034 *down_cast<Geographic_multilinestring *>(g1_mls.get()), *g2,
1035 m_geographic_ll_la_aa_strategy)) &&
1036 (g1_mpy->empty() ||
1037 bg::covered_by(*down_cast<Geographic_multipolygon *>(g1_mpy.get()),
1038 *g2, m_geographic_ll_la_aa_strategy));
1039 }
1040 if (eval(down_cast<Geographic_multilinestring *>(g1_mls.get()), g2)) {
1041 for (auto &pt : *down_cast<Geographic_multipoint *>(g1_mpt.get()))
1042 if (!bg::covered_by(pt, *g2, m_geographic_pl_pa_strategy)) return false;
1043 return g1_mpy->empty() ||
1044 bg::covered_by(*down_cast<Geographic_multipolygon *>(g1_mpy.get()),
1045 *g2, m_geographic_ll_la_aa_strategy);
1046 }
1047 if (eval(down_cast<Geographic_multipolygon *>(g1_mpy.get()), g2)) {
1048 for (auto &pt : *down_cast<Geographic_multipoint *>(g1_mpt.get()))
1049 if (!bg::covered_by(pt, *g2, m_geographic_pl_pa_strategy)) return false;
1050 return g1_mls->empty() ||
1051 bg::covered_by(
1052 *down_cast<Geographic_multilinestring *>(g1_mls.get()), *g2,
1053 m_geographic_ll_la_aa_strategy);
1054 }
1055 return false;
1056 }
1057
1058 //////////////////////////////////////////////////////////////////////////////
1059
1060 // within(Geographic_multipoint, *)
1061
eval(const Geographic_multipoint * g1,const Geographic_point * g2) const1062 bool Within::eval(const Geographic_multipoint *g1,
1063 const Geographic_point *g2) const {
1064 Equals equals(m_semi_major, m_semi_minor);
1065 return equals(g1, g2);
1066 }
1067
eval(const Geographic_multipoint * g1,const Geographic_linestring * g2) const1068 bool Within::eval(const Geographic_multipoint *g1,
1069 const Geographic_linestring *g2) const {
1070 // At least one point in g1 must be within g2. The rest has to intersect g2.
1071 bool within = false;
1072 bool intersects = false;
1073 for (auto &pt : *g1) {
1074 if (!within) {
1075 within = bg::within(pt, *g2, m_geographic_pl_pa_strategy);
1076 if (!within)
1077 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1078 else
1079 intersects = true;
1080 } else {
1081 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1082 }
1083 if (!intersects) break;
1084 }
1085 return (within && intersects);
1086 }
1087
eval(const Geographic_multipoint * g1,const Geographic_polygon * g2) const1088 bool Within::eval(const Geographic_multipoint *g1,
1089 const Geographic_polygon *g2) const {
1090 // At least one point in g1 must be within g2. The rest has to intersect g2.
1091 bool within = false;
1092 bool intersects = false;
1093 for (auto &pt : *g1) {
1094 if (!within) {
1095 within = bg::within(pt, *g2, m_geographic_pl_pa_strategy);
1096 if (!within)
1097 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1098 else
1099 intersects = true;
1100 } else {
1101 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1102 }
1103 if (!intersects) break;
1104 }
1105 return (within && intersects);
1106 }
1107
eval(const Geographic_multipoint * g1,const Geographic_geometrycollection * g2) const1108 bool Within::eval(const Geographic_multipoint *g1,
1109 const Geographic_geometrycollection *g2) const {
1110 // At least one point in g1 must be within g2. The rest has to intersect g2.
1111 Intersects intersects_func(m_semi_major, m_semi_minor);
1112 bool within = false;
1113 bool intersects = false;
1114 for (auto &pt : *g1) {
1115 if (!within) {
1116 within = eval(&pt, g2);
1117 if (!within)
1118 intersects = intersects_func(&pt, g2);
1119 else
1120 intersects = true;
1121 } else {
1122 intersects = intersects_func(&pt, g2);
1123 }
1124 if (!intersects) break;
1125 }
1126 return (within && intersects);
1127 }
1128
eval(const Geographic_multipoint * g1,const Geographic_multipoint * g2) const1129 bool Within::eval(const Geographic_multipoint *g1,
1130 const Geographic_multipoint *g2) const {
1131 // Default strategy is OK. P/P computations do not depend on shape of
1132 // ellipsoid.
1133 for (auto &pt : *g1)
1134 if (!bg::within(pt, *g2)) return false;
1135 return true;
1136 }
1137
eval(const Geographic_multipoint * g1,const Geographic_multilinestring * g2) const1138 bool Within::eval(const Geographic_multipoint *g1,
1139 const Geographic_multilinestring *g2) const {
1140 // At least one point in g1 must be within g2. The rest has to intersect g2.
1141 bool within = false;
1142 bool intersects = false;
1143 for (auto &pt : *g1) {
1144 if (!within) {
1145 within = bg::within(pt, *g2, m_geographic_pl_pa_strategy);
1146 if (!within)
1147 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1148 else
1149 intersects = true;
1150 } else {
1151 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1152 }
1153 if (!intersects) break;
1154 }
1155 return (within && intersects);
1156 }
1157
eval(const Geographic_multipoint * g1,const Geographic_multipolygon * g2) const1158 bool Within::eval(const Geographic_multipoint *g1,
1159 const Geographic_multipolygon *g2) const {
1160 // At least one point in g1 must be within g2. The rest has to intersect g2.
1161 bool within = false;
1162 bool intersects = false;
1163 for (auto &pt : *g1) {
1164 if (!within) {
1165 within = bg::within(pt, *g2, m_geographic_pl_pa_strategy);
1166 if (!within)
1167 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1168 else
1169 intersects = true;
1170 } else {
1171 intersects = bg::intersects(pt, *g2, m_geographic_pl_pa_strategy);
1172 }
1173 if (!intersects) break;
1174 }
1175 return (within && intersects);
1176 }
1177
1178 //////////////////////////////////////////////////////////////////////////////
1179
1180 // within(Geographic_multilinestring, *)
1181
eval(const Geographic_multilinestring *,const Geographic_point *) const1182 bool Within::eval(const Geographic_multilinestring *,
1183 const Geographic_point *) const {
1184 // A multilinestring can never be within a point.
1185 return false;
1186 }
1187
eval(const Geographic_multilinestring * g1,const Geographic_linestring * g2) const1188 bool Within::eval(const Geographic_multilinestring *g1,
1189 const Geographic_linestring *g2) const {
1190 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
1191 }
1192
eval(const Geographic_multilinestring * g1,const Geographic_polygon * g2) const1193 bool Within::eval(const Geographic_multilinestring *g1,
1194 const Geographic_polygon *g2) const {
1195 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
1196 }
1197
eval(const Geographic_multilinestring * g1,const Geographic_geometrycollection * g2) const1198 bool Within::eval(const Geographic_multilinestring *g1,
1199 const Geographic_geometrycollection *g2) const {
1200 // For g1 to be within g2, no point of g1 may be in the exterior of g2 and at
1201 // least one point of the interior of g1 has to be within the interior of g2.
1202
1203 std::unique_ptr<Multipoint> g2_mpt;
1204 std::unique_ptr<Multilinestring> g2_mls;
1205 std::unique_ptr<Multipolygon> g2_mpy;
1206 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
1207 &g2_mpy);
1208 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
1209
1210 Difference difference(m_semi_major, m_semi_minor);
1211 std::unique_ptr<Geographic_multilinestring> g1_diff_g2(
1212 new Geographic_multilinestring());
1213 g1_diff_g2.reset(down_cast<Geographic_multilinestring *>(
1214 difference(g1, down_cast<Geographic_multilinestring *>(g2_mls.get()))));
1215 g1_diff_g2.reset(down_cast<Geographic_multilinestring *>(difference(
1216 g1_diff_g2.get(), down_cast<Geographic_multipolygon *>(g2_mpy.get()))));
1217
1218 boost::geometry::de9im::mask mask("T********");
1219 return g1_diff_g2->empty() &&
1220 (bg::relate(*g1,
1221 *down_cast<Geographic_multilinestring *>(g2_mls.get()),
1222 mask, m_geographic_ll_la_aa_strategy) ||
1223 bg::relate(*g1, *down_cast<Geographic_multipolygon *>(g2_mpy.get()),
1224 mask, m_geographic_ll_la_aa_strategy));
1225 }
1226
eval(const Geographic_multilinestring *,const Geographic_multipoint *) const1227 bool Within::eval(const Geographic_multilinestring *,
1228 const Geographic_multipoint *) const {
1229 // A multilinestring can never be within a multipoint.
1230 return false;
1231 }
1232
eval(const Geographic_multilinestring * g1,const Geographic_multilinestring * g2) const1233 bool Within::eval(const Geographic_multilinestring *g1,
1234 const Geographic_multilinestring *g2) const {
1235 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
1236 }
1237
eval(const Geographic_multilinestring * g1,const Geographic_multipolygon * g2) const1238 bool Within::eval(const Geographic_multilinestring *g1,
1239 const Geographic_multipolygon *g2) const {
1240 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
1241 }
1242
1243 //////////////////////////////////////////////////////////////////////////////
1244
1245 // within(Geographic_multipolygon, *)
1246
eval(const Geographic_multipolygon *,const Geographic_point *) const1247 bool Within::eval(const Geographic_multipolygon *,
1248 const Geographic_point *) const {
1249 // A multipolygon can never be within a point.
1250 return false;
1251 }
1252
eval(const Geographic_multipolygon *,const Geographic_linestring *) const1253 bool Within::eval(const Geographic_multipolygon *,
1254 const Geographic_linestring *) const {
1255 // A multipolygon can never be within a linestring.
1256 return false;
1257 }
1258
eval(const Geographic_multipolygon * g1,const Geographic_polygon * g2) const1259 bool Within::eval(const Geographic_multipolygon *g1,
1260 const Geographic_polygon *g2) const {
1261 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
1262 }
1263
eval(const Geographic_multipolygon * g1,const Geographic_geometrycollection * g2) const1264 bool Within::eval(const Geographic_multipolygon *g1,
1265 const Geographic_geometrycollection *g2) const {
1266 // A multipolygon may not be within the points and linestrings of g2, so the
1267 // only way a multipolygon is within a geometrycollectin, is if it's within
1268 // the union multipolygon of the geometrycollection.
1269 std::unique_ptr<Multipoint> g2_mpt;
1270 std::unique_ptr<Multilinestring> g2_mls;
1271 std::unique_ptr<Multipolygon> g2_mpy;
1272 split_gc(down_cast<const Geometrycollection *>(g2), &g2_mpt, &g2_mls,
1273 &g2_mpy);
1274 gc_union(m_semi_major, m_semi_minor, &g2_mpt, &g2_mls, &g2_mpy);
1275 return eval(g1, down_cast<Geographic_multipolygon *>(g2_mpy.get()));
1276 }
1277
eval(const Geographic_multipolygon *,const Geographic_multipoint *) const1278 bool Within::eval(const Geographic_multipolygon *,
1279 const Geographic_multipoint *) const {
1280 // A multipolygon can never be within a multipoint.
1281 return false;
1282 }
1283
eval(const Geographic_multipolygon *,const Geographic_multilinestring *) const1284 bool Within::eval(const Geographic_multipolygon *,
1285 const Geographic_multilinestring *) const {
1286 // A multipolygon can never be within a multilinestring.
1287 return false;
1288 }
1289
eval(const Geographic_multipolygon * g1,const Geographic_multipolygon * g2) const1290 bool Within::eval(const Geographic_multipolygon *g1,
1291 const Geographic_multipolygon *g2) const {
1292 return bg::within(*g1, *g2, m_geographic_ll_la_aa_strategy);
1293 }
1294
1295 //////////////////////////////////////////////////////////////////////////////
1296
1297 // within(Box, Box)
1298
eval(const Cartesian_box * b1,const Cartesian_box * b2) const1299 bool Within::eval(const Cartesian_box *b1, const Cartesian_box *b2) const {
1300 if (mbrs_are_equal(*b1, *b2)) return true;
1301
1302 // Work around bugs in BG for boxes that have zero height and/or width.
1303 if (mbr_is_point(*b1)) {
1304 Cartesian_point pt(b1->min_corner().x(), b1->min_corner().y());
1305
1306 if (mbr_is_line(*b2)) {
1307 Cartesian_point b2_ls_start(b2->min_corner().x(), b2->min_corner().y());
1308 Cartesian_point b2_ls_end(b2->max_corner().x(), b2->max_corner().y());
1309 Cartesian_linestring b2_ls;
1310 b2_ls.push_back(b2_ls_start);
1311 b2_ls.push_back(b2_ls_end);
1312
1313 return bg::within(pt, b2_ls);
1314 }
1315
1316 return bg::within(pt, *b2);
1317 }
1318
1319 if (mbr_is_line(*b1)) {
1320 Cartesian_point b1_ls_start(b1->min_corner().x(), b1->min_corner().y());
1321 Cartesian_point b1_ls_end(b1->max_corner().x(), b1->max_corner().y());
1322 Cartesian_linestring b1_ls;
1323 b1_ls.push_back(b1_ls_start);
1324 b1_ls.push_back(b1_ls_end);
1325
1326 if (mbr_is_line(*b2)) {
1327 Cartesian_point b2_ls_start(b2->min_corner().x(), b2->min_corner().y());
1328 Cartesian_point b2_ls_end(b2->max_corner().x(), b2->max_corner().y());
1329 Cartesian_linestring b2_ls;
1330 b2_ls.push_back(b2_ls_start);
1331 b2_ls.push_back(b2_ls_end);
1332
1333 return bg::within(b1_ls, b2_ls);
1334 }
1335
1336 Cartesian_point b2_pt1(b2->min_corner().x(), b2->min_corner().y());
1337 Cartesian_point b2_pt2(b2->max_corner().x(), b2->min_corner().y());
1338 Cartesian_point b2_pt3(b2->max_corner().x(), b2->max_corner().y());
1339 Cartesian_point b2_pt4(b2->min_corner().x(), b2->max_corner().y());
1340 Cartesian_point b2_pt5(b2->min_corner().x(), b2->min_corner().y());
1341 Cartesian_linearring b2_lr;
1342 b2_lr.push_back(b2_pt1);
1343 b2_lr.push_back(b2_pt2);
1344 b2_lr.push_back(b2_pt3);
1345 b2_lr.push_back(b2_pt4);
1346 b2_lr.push_back(b2_pt5);
1347 Cartesian_polygon b2_py;
1348 b2_py.push_back(b2_lr);
1349
1350 return bg::within(b1_ls, b2_py);
1351 }
1352
1353 return bg::within(*b1, *b2);
1354 }
1355
eval(const Geographic_box * b1,const Geographic_box * b2) const1356 bool Within::eval(const Geographic_box *b1, const Geographic_box *b2) const {
1357 if (mbrs_are_equal(*b1, *b2)) return true;
1358
1359 // Work around bugs in BG for boxes that have zero height and/or width.
1360 if (mbr_is_point(*b1)) {
1361 Geographic_point pt(b1->min_corner().x(), b1->min_corner().y());
1362
1363 if (mbr_is_line(*b2)) {
1364 Geographic_point b2_ls_start(b2->min_corner().x(), b2->min_corner().y());
1365 Geographic_point b2_ls_end(b2->max_corner().x(), b2->max_corner().y());
1366 Geographic_linestring b2_ls;
1367 b2_ls.push_back(b2_ls_start);
1368 b2_ls.push_back(b2_ls_end);
1369
1370 return bg::within(pt, b2_ls);
1371 }
1372
1373 return bg::within(pt, *b2);
1374 }
1375
1376 if (mbr_is_line(*b1)) {
1377 Geographic_point b1_ls_start(b1->min_corner().x(), b1->min_corner().y());
1378 Geographic_point b1_ls_end(b1->max_corner().x(), b1->max_corner().y());
1379 Geographic_linestring b1_ls;
1380 b1_ls.push_back(b1_ls_start);
1381 b1_ls.push_back(b1_ls_end);
1382
1383 if (mbr_is_line(*b2)) {
1384 Geographic_point b2_ls_start(b2->min_corner().x(), b2->min_corner().y());
1385 Geographic_point b2_ls_end(b2->max_corner().x(), b2->max_corner().y());
1386 Geographic_linestring b2_ls;
1387 b2_ls.push_back(b2_ls_start);
1388 b2_ls.push_back(b2_ls_end);
1389
1390 return bg::within(b1_ls, b2_ls);
1391 }
1392
1393 // If b1 is a line along the border of b2, it's not within b2.
1394 if (((b1_ls_start.x() == b1_ls_end.x()) &&
1395 (b1_ls_start.x() == b2->min_corner().x() ||
1396 b1_ls_start.x() == b2->max_corner().x())) ||
1397 ((b1_ls_start.y() == b1_ls_end.y()) &&
1398 (b1_ls_start.y() == b2->min_corner().y() ||
1399 b1_ls_start.y() == b2->max_corner().y())))
1400 return false;
1401
1402 return bg::covered_by(b1_ls_start, *b2) && bg::covered_by(b1_ls_end, *b2);
1403 }
1404
1405 return bg::within(*b1, *b2);
1406 }
1407
1408 //////////////////////////////////////////////////////////////////////////////
1409
within(const dd::Spatial_reference_system * srs,const Geometry * g1,const Geometry * g2,const char * func_name,bool * within,bool * null)1410 bool within(const dd::Spatial_reference_system *srs, const Geometry *g1,
1411 const Geometry *g2, const char *func_name, bool *within,
1412 bool *null) noexcept {
1413 try {
1414 DBUG_ASSERT(g1->coordinate_system() == g2->coordinate_system());
1415 DBUG_ASSERT(srs == nullptr ||
1416 ((srs->is_cartesian() &&
1417 g1->coordinate_system() == Coordinate_system::kCartesian) ||
1418 (srs->is_geographic() &&
1419 g1->coordinate_system() == Coordinate_system::kGeographic)));
1420
1421 if ((*null = (g1->is_empty() || g2->is_empty()))) return false;
1422
1423 Within within_func(srs ? srs->semi_major_axis() : 0.0,
1424 srs ? srs->semi_minor_axis() : 0.0);
1425 *within = within_func(g1, g2);
1426 } catch (...) {
1427 handle_gis_exception(func_name);
1428 return true;
1429 }
1430
1431 return false;
1432 }
1433
mbr_within(const dd::Spatial_reference_system * srs,const Geometry * g1,const Geometry * g2,const char * func_name,bool * within,bool * null)1434 bool mbr_within(const dd::Spatial_reference_system *srs, const Geometry *g1,
1435 const Geometry *g2, const char *func_name, bool *within,
1436 bool *null) noexcept {
1437 try {
1438 DBUG_ASSERT(g1->coordinate_system() == g2->coordinate_system());
1439 DBUG_ASSERT(srs == nullptr ||
1440 ((srs->is_cartesian() &&
1441 g1->coordinate_system() == Coordinate_system::kCartesian) ||
1442 (srs->is_geographic() &&
1443 g1->coordinate_system() == Coordinate_system::kGeographic)));
1444
1445 if ((*null = (g1->is_empty() || g2->is_empty()))) return false;
1446
1447 Within within_func(srs ? srs->semi_major_axis() : 0.0,
1448 srs ? srs->semi_minor_axis() : 0.0);
1449
1450 switch (g1->coordinate_system()) {
1451 case Coordinate_system::kCartesian: {
1452 Cartesian_box mbr1;
1453 box_envelope(g1, srs, &mbr1);
1454 Cartesian_box mbr2;
1455 box_envelope(g2, srs, &mbr2);
1456
1457 *within = within_func(&mbr1, &mbr2);
1458 break;
1459 }
1460 case Coordinate_system::kGeographic: {
1461 Geographic_box mbr1;
1462 box_envelope(g1, srs, &mbr1);
1463 Geographic_box mbr2;
1464 box_envelope(g2, srs, &mbr2);
1465
1466 *within = within_func(&mbr1, &mbr2);
1467 break;
1468 }
1469 }
1470 } catch (...) {
1471 handle_gis_exception(func_name);
1472 return true;
1473 }
1474
1475 return false;
1476 }
1477
1478 } // namespace gis
1479