1 | #pragma once |
2 | |
3 | #include <Core/Types.h> |
4 | #include <Core/Defines.h> |
5 | #include <Core/TypeListNumber.h> |
6 | #include <Columns/IColumn.h> |
7 | #include <Columns/ColumnVector.h> |
8 | #include <Common/typeid_cast.h> |
9 | #include <ext/range.h> |
10 | |
11 | /// Warning in boost::geometry during template strategy substitution. |
12 | #pragma GCC diagnostic push |
13 | |
14 | #if !__clang__ |
15 | #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" |
16 | #endif |
17 | |
18 | #pragma GCC diagnostic ignored "-Wunused-parameter" |
19 | |
20 | #include <boost/geometry.hpp> |
21 | |
22 | #pragma GCC diagnostic pop |
23 | |
24 | #include <boost/geometry/geometries/point_xy.hpp> |
25 | #include <boost/geometry/geometries/polygon.hpp> |
26 | #include <boost/geometry/geometries/multi_polygon.hpp> |
27 | #include <boost/geometry/geometries/segment.hpp> |
28 | #include <boost/geometry/algorithms/comparable_distance.hpp> |
29 | #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp> |
30 | |
31 | #include <array> |
32 | #include <vector> |
33 | #include <iterator> |
34 | #include <cmath> |
35 | #include <algorithm> |
36 | #include <IO/WriteBufferFromString.h> |
37 | |
38 | namespace DB |
39 | { |
40 | |
41 | namespace ErrorCodes |
42 | { |
43 | extern const int LOGICAL_ERROR; |
44 | } |
45 | |
46 | |
47 | namespace GeoUtils |
48 | { |
49 | |
50 | |
51 | template <typename Polygon> |
52 | UInt64 getPolygonAllocatedBytes(const Polygon & polygon) |
53 | { |
54 | UInt64 size = 0; |
55 | |
56 | using RingType = typename Polygon::ring_type; |
57 | using ValueType = typename RingType::value_type; |
58 | |
59 | auto sizeOfRing = [](const RingType & ring) { return sizeof(ring) + ring.capacity() * sizeof(ValueType); }; |
60 | |
61 | size += sizeOfRing(polygon.outer()); |
62 | |
63 | const auto & inners = polygon.inners(); |
64 | size += sizeof(inners) + inners.capacity() * sizeof(RingType); |
65 | for (auto & inner : inners) |
66 | size += sizeOfRing(inner); |
67 | |
68 | return size; |
69 | } |
70 | |
71 | template <typename MultiPolygon> |
72 | UInt64 getMultiPolygonAllocatedBytes(const MultiPolygon & multi_polygon) |
73 | { |
74 | using ValueType = typename MultiPolygon::value_type; |
75 | UInt64 size = multi_polygon.capacity() * sizeof(ValueType); |
76 | |
77 | for (const auto & polygon : multi_polygon) |
78 | size += getPolygonAllocatedBytes(polygon); |
79 | |
80 | return size; |
81 | } |
82 | |
83 | template <typename CoordinateType = Float32> |
84 | class PointInPolygonWithGrid |
85 | { |
86 | public: |
87 | using Point = boost::geometry::model::d2::point_xy<CoordinateType>; |
88 | /// Counter-Clockwise ordering. |
89 | using Polygon = boost::geometry::model::polygon<Point, false>; |
90 | using MultiPolygon = boost::geometry::model::multi_polygon<Polygon>; |
91 | using Box = boost::geometry::model::box<Point>; |
92 | using Segment = boost::geometry::model::segment<Point>; |
93 | |
94 | explicit PointInPolygonWithGrid(const Polygon & polygon_, UInt16 grid_size_ = 8) |
95 | : grid_size(std::max<UInt16>(1, grid_size_)), polygon(polygon_) {} |
96 | |
97 | void init(); |
98 | |
99 | /// True if bound box is empty. |
100 | bool hasEmptyBound() const { return has_empty_bound; } |
101 | |
102 | UInt64 getAllocatedBytes() const; |
103 | |
104 | inline bool ALWAYS_INLINE contains(CoordinateType x, CoordinateType y); |
105 | |
106 | private: |
107 | enum class CellType |
108 | { |
109 | inner, |
110 | outer, |
111 | singleLine, |
112 | pairOfLinesSingleConvexPolygon, |
113 | pairOfLinesSingleNonConvexPolygons, |
114 | pairOfLinesDifferentPolygons, |
115 | complexPolygon |
116 | }; |
117 | |
118 | struct HalfPlane |
119 | { |
120 | /// Line, a * x + b * y + c = 0. Vector (a, b) points inside half-plane. |
121 | CoordinateType a; |
122 | CoordinateType b; |
123 | CoordinateType c; |
124 | |
125 | /// Take left half-plane. |
126 | void fill(const Point & from, const Point & to) |
127 | { |
128 | a = -(to.y() - from.y()); |
129 | b = to.x() - from.x(); |
130 | c = -from.x() * a - from.y() * b; |
131 | } |
132 | |
133 | /// Inner part of the HalfPlane is the left side of initialized vector. |
134 | bool ALWAYS_INLINE contains(CoordinateType x, CoordinateType y) const { return a * x + b * y + c >= 0; } |
135 | }; |
136 | |
137 | struct Cell |
138 | { |
139 | static const int max_stored_half_planes = 2; |
140 | |
141 | HalfPlane half_planes[max_stored_half_planes]; |
142 | size_t index_of_inner_polygon; |
143 | CellType type; |
144 | }; |
145 | |
146 | const UInt16 grid_size; |
147 | |
148 | Polygon polygon; |
149 | std::vector<Cell> cells; |
150 | std::vector<MultiPolygon> polygons; |
151 | |
152 | CoordinateType cell_width; |
153 | CoordinateType cell_height; |
154 | |
155 | CoordinateType x_shift; |
156 | CoordinateType y_shift; |
157 | CoordinateType x_scale; |
158 | CoordinateType y_scale; |
159 | |
160 | bool has_empty_bound = false; |
161 | bool was_grid_built = false; |
162 | |
163 | void buildGrid(); |
164 | |
165 | void calcGridAttributes(Box & box); |
166 | |
167 | template <typename T> |
168 | T ALWAYS_INLINE getCellIndex(T row, T col) const { return row * grid_size + col; } |
169 | |
170 | /// Complex case. Will check intersection directly. |
171 | inline void addComplexPolygonCell(size_t index, const Box & box); |
172 | |
173 | /// Empty intersection or intersection == box. |
174 | inline void addCell(size_t index, const Box & empty_box); |
175 | |
176 | /// Intersection is a single polygon. |
177 | inline void addCell(size_t index, const Box & box, const Polygon & intersection); |
178 | |
179 | /// Intersection is a pair of polygons. |
180 | inline void addCell(size_t index, const Box & box, const Polygon & first, const Polygon & second); |
181 | |
182 | /// Returns a list of half-planes were formed from intersection edges without box edges. |
183 | inline std::vector<HalfPlane> findHalfPlanes(const Box & box, const Polygon & intersection); |
184 | |
185 | /// Check that polygon.outer() is convex. |
186 | inline bool isConvex(const Polygon & polygon); |
187 | |
188 | using Distance = typename boost::geometry::default_comparable_distance_result<Point, Segment>::type; |
189 | |
190 | /// min(distance(point, edge) : edge in polygon) |
191 | inline Distance distance(const Point & point, const Polygon & polygon); |
192 | }; |
193 | |
194 | template <typename CoordinateType> |
195 | UInt64 PointInPolygonWithGrid<CoordinateType>::getAllocatedBytes() const |
196 | { |
197 | UInt64 size = sizeof(*this); |
198 | |
199 | size += cells.capacity() * sizeof(Cell); |
200 | size += polygons.capacity() * sizeof(MultiPolygon); |
201 | size += getPolygonAllocatedBytes(polygon); |
202 | |
203 | for (const auto & elem : polygons) |
204 | size += getMultiPolygonAllocatedBytes(elem); |
205 | |
206 | return size; |
207 | } |
208 | |
209 | template <typename CoordinateType> |
210 | void PointInPolygonWithGrid<CoordinateType>::init() |
211 | { |
212 | if (!was_grid_built) |
213 | buildGrid(); |
214 | |
215 | was_grid_built = true; |
216 | } |
217 | |
218 | template <typename CoordinateType> |
219 | void PointInPolygonWithGrid<CoordinateType>::calcGridAttributes( |
220 | PointInPolygonWithGrid<CoordinateType>::Box & box) |
221 | { |
222 | boost::geometry::envelope(polygon, box); |
223 | |
224 | const Point & min_corner = box.min_corner(); |
225 | const Point & max_corner = box.max_corner(); |
226 | |
227 | #pragma GCC diagnostic push |
228 | #if !__clang__ |
229 | #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" |
230 | #endif |
231 | |
232 | cell_width = (max_corner.x() - min_corner.x()) / grid_size; |
233 | cell_height = (max_corner.y() - min_corner.y()) / grid_size; |
234 | |
235 | #pragma GCC diagnostic pop |
236 | |
237 | if (cell_width == 0 || cell_height == 0) |
238 | { |
239 | has_empty_bound = true; |
240 | return; |
241 | } |
242 | |
243 | x_scale = 1 / cell_width; |
244 | y_scale = 1 / cell_height; |
245 | x_shift = -min_corner.x(); |
246 | y_shift = -min_corner.y(); |
247 | } |
248 | |
249 | template <typename CoordinateType> |
250 | void PointInPolygonWithGrid<CoordinateType>::buildGrid() |
251 | { |
252 | Box box; |
253 | calcGridAttributes(box); |
254 | |
255 | if (has_empty_bound) |
256 | return; |
257 | |
258 | cells.assign(grid_size * grid_size, {}); |
259 | |
260 | const Point & min_corner = box.min_corner(); |
261 | |
262 | for (size_t row = 0; row < grid_size; ++row) |
263 | { |
264 | #pragma GCC diagnostic push |
265 | #if !__clang__ |
266 | #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" |
267 | #endif |
268 | CoordinateType y_min = min_corner.y() + row * cell_height; |
269 | CoordinateType y_max = min_corner.y() + (row + 1) * cell_height; |
270 | |
271 | for (size_t col = 0; col < grid_size; ++col) |
272 | { |
273 | CoordinateType x_min = min_corner.x() + col * cell_width; |
274 | CoordinateType x_max = min_corner.x() + (col + 1) * cell_width; |
275 | #pragma GCC diagnostic pop |
276 | Box cell_box(Point(x_min, y_min), Point(x_max, y_max)); |
277 | |
278 | Polygon cell_bound; |
279 | boost::geometry::convert(cell_box, cell_bound); |
280 | |
281 | MultiPolygon intersection; |
282 | boost::geometry::intersection(polygon, cell_bound, intersection); |
283 | |
284 | size_t cellIndex = getCellIndex(row, col); |
285 | |
286 | if (intersection.empty()) |
287 | addCell(cellIndex, cell_box); |
288 | else if (intersection.size() == 1) |
289 | addCell(cellIndex, cell_box, intersection.front()); |
290 | else if (intersection.size() == 2) |
291 | addCell(cellIndex, cell_box, intersection.front(), intersection.back()); |
292 | else |
293 | addComplexPolygonCell(cellIndex, cell_box); |
294 | } |
295 | } |
296 | } |
297 | |
298 | template <typename CoordinateType> |
299 | bool PointInPolygonWithGrid<CoordinateType>::contains(CoordinateType x, CoordinateType y) |
300 | { |
301 | if (has_empty_bound) |
302 | return false; |
303 | |
304 | CoordinateType float_row = (y + y_shift) * y_scale; |
305 | CoordinateType float_col = (x + x_shift) * x_scale; |
306 | |
307 | if (float_row < 0 || float_row > grid_size) |
308 | return false; |
309 | if (float_col < 0 || float_col > grid_size) |
310 | return false; |
311 | |
312 | int row = std::min<int>(float_row, grid_size - 1); |
313 | int col = std::min<int>(float_col, grid_size - 1); |
314 | |
315 | int index = getCellIndex(row, col); |
316 | const auto & cell = cells[index]; |
317 | |
318 | switch (cell.type) |
319 | { |
320 | case CellType::inner: |
321 | return true; |
322 | case CellType::outer: |
323 | return false; |
324 | case CellType::singleLine: |
325 | return cell.half_planes[0].contains(x, y); |
326 | case CellType::pairOfLinesSingleConvexPolygon: |
327 | return cell.half_planes[0].contains(x, y) && cell.half_planes[1].contains(x, y); |
328 | case CellType::pairOfLinesDifferentPolygons: [[fallthrough]]; |
329 | case CellType::pairOfLinesSingleNonConvexPolygons: |
330 | return cell.half_planes[0].contains(x, y) || cell.half_planes[1].contains(x, y); |
331 | case CellType::complexPolygon: |
332 | return boost::geometry::within(Point(x, y), polygons[cell.index_of_inner_polygon]); |
333 | } |
334 | |
335 | __builtin_unreachable(); |
336 | } |
337 | |
338 | template <typename CoordinateType> |
339 | typename PointInPolygonWithGrid<CoordinateType>::Distance |
340 | PointInPolygonWithGrid<CoordinateType>::distance( |
341 | const PointInPolygonWithGrid<CoordinateType>::Point & point, |
342 | const PointInPolygonWithGrid<CoordinateType>::Polygon & poly) |
343 | { |
344 | const auto & outer = poly.outer(); |
345 | Distance distance = 0; |
346 | for (auto i : ext::range(0, outer.size() - 1)) |
347 | { |
348 | Segment segment(outer[i], outer[i + 1]); |
349 | Distance current = boost::geometry::comparable_distance(point, segment); |
350 | distance = i ? std::min(current, distance) : current; |
351 | } |
352 | return distance; |
353 | } |
354 | |
355 | template <typename CoordinateType> |
356 | bool PointInPolygonWithGrid<CoordinateType>::isConvex(const PointInPolygonWithGrid<CoordinateType>::Polygon & poly) |
357 | { |
358 | const auto & outer = poly.outer(); |
359 | /// Segment or point. |
360 | if (outer.size() < 4) |
361 | return false; |
362 | |
363 | auto vecProduct = [](const Point & from, const Point & to) { return from.x() * to.y() - from.y() * to.x(); }; |
364 | auto getVector = [](const Point & from, const Point & to) -> Point |
365 | { |
366 | return Point(to.x() - from.x(), to.y() - from.y()); |
367 | }; |
368 | |
369 | Point first = getVector(outer[0], outer[1]); |
370 | Point prev = first; |
371 | |
372 | for (auto i : ext::range(1, outer.size() - 1)) |
373 | { |
374 | Point cur = getVector(outer[i], outer[i + 1]); |
375 | if (vecProduct(prev, cur) < 0) |
376 | return false; |
377 | |
378 | prev = cur; |
379 | } |
380 | |
381 | return vecProduct(prev, first) >= 0; |
382 | } |
383 | |
384 | template <typename CoordinateType> |
385 | std::vector<typename PointInPolygonWithGrid<CoordinateType>::HalfPlane> |
386 | PointInPolygonWithGrid<CoordinateType>::findHalfPlanes( |
387 | const PointInPolygonWithGrid<CoordinateType>::Box & box, |
388 | const PointInPolygonWithGrid<CoordinateType>::Polygon & intersection) |
389 | { |
390 | std::vector<HalfPlane> half_planes; |
391 | Polygon bound; |
392 | boost::geometry::convert(box, bound); |
393 | const auto & outer = intersection.outer(); |
394 | |
395 | for (auto i : ext::range(0, outer.size() - 1)) |
396 | { |
397 | /// Want to detect is intersection edge was formed from box edge or from polygon edge. |
398 | /// If center of the edge closer to box, than don't form the half-plane. |
399 | Segment segment(outer[i], outer[i + 1]); |
400 | Point center((segment.first.x() + segment.second.x()) / 2, (segment.first.y() + segment.second.y()) / 2); |
401 | if (distance(center, polygon) < distance(center, bound)) |
402 | { |
403 | half_planes.push_back({}); |
404 | half_planes.back().fill(segment.first, segment.second); |
405 | } |
406 | } |
407 | |
408 | return half_planes; |
409 | } |
410 | |
411 | template <typename CoordinateType> |
412 | void PointInPolygonWithGrid<CoordinateType>::addComplexPolygonCell( |
413 | size_t index, const PointInPolygonWithGrid<CoordinateType>::Box & box) |
414 | { |
415 | cells[index].type = CellType::complexPolygon; |
416 | cells[index].index_of_inner_polygon = polygons.size(); |
417 | |
418 | /// Expand box in (1 + eps_factor) times to eliminate errors for points on box bound. |
419 | static constexpr CoordinateType eps_factor = 0.01; |
420 | auto x_eps = eps_factor * (box.max_corner().x() - box.min_corner().x()); |
421 | auto y_eps = eps_factor * (box.max_corner().y() - box.min_corner().y()); |
422 | |
423 | Point min_corner(box.min_corner().x() - x_eps, box.min_corner().y() - y_eps); |
424 | Point max_corner(box.max_corner().x() + x_eps, box.max_corner().y() + y_eps); |
425 | Box box_with_eps_bound(min_corner, max_corner); |
426 | |
427 | Polygon bound; |
428 | boost::geometry::convert(box_with_eps_bound, bound); |
429 | |
430 | MultiPolygon intersection; |
431 | boost::geometry::intersection(polygon, bound, intersection); |
432 | |
433 | polygons.push_back(intersection); |
434 | } |
435 | |
436 | template <typename CoordinateType> |
437 | void PointInPolygonWithGrid<CoordinateType>::addCell( |
438 | size_t index, const PointInPolygonWithGrid<CoordinateType>::Box & empty_box) |
439 | { |
440 | const auto & min_corner = empty_box.min_corner(); |
441 | const auto & max_corner = empty_box.max_corner(); |
442 | |
443 | Point center((min_corner.x() + max_corner.x()) / 2, (min_corner.y() + max_corner.y()) / 2); |
444 | |
445 | if (boost::geometry::within(center, polygon)) |
446 | cells[index].type = CellType::inner; |
447 | else |
448 | cells[index].type = CellType::outer; |
449 | |
450 | } |
451 | |
452 | template <typename CoordinateType> |
453 | void PointInPolygonWithGrid<CoordinateType>::addCell( |
454 | size_t index, |
455 | const PointInPolygonWithGrid<CoordinateType>::Box & box, |
456 | const PointInPolygonWithGrid<CoordinateType>::Polygon & intersection) |
457 | { |
458 | if (!intersection.inners().empty()) |
459 | addComplexPolygonCell(index, box); |
460 | |
461 | auto half_planes = findHalfPlanes(box, intersection); |
462 | |
463 | if (half_planes.empty()) |
464 | addCell(index, box); |
465 | else if (half_planes.size() == 1) |
466 | { |
467 | cells[index].type = CellType::singleLine; |
468 | cells[index].half_planes[0] = half_planes[0]; |
469 | } |
470 | else if (half_planes.size() == 2) |
471 | { |
472 | cells[index].type = isConvex(intersection) ? CellType::pairOfLinesSingleConvexPolygon |
473 | : CellType::pairOfLinesSingleNonConvexPolygons; |
474 | cells[index].half_planes[0] = half_planes[0]; |
475 | cells[index].half_planes[1] = half_planes[1]; |
476 | } |
477 | else |
478 | addComplexPolygonCell(index, box); |
479 | } |
480 | |
481 | template <typename CoordinateType> |
482 | void PointInPolygonWithGrid<CoordinateType>::addCell( |
483 | size_t index, |
484 | const PointInPolygonWithGrid<CoordinateType>::Box & box, |
485 | const PointInPolygonWithGrid<CoordinateType>::Polygon & first, |
486 | const PointInPolygonWithGrid<CoordinateType>::Polygon & second) |
487 | { |
488 | if (!first.inners().empty() || !second.inners().empty()) |
489 | addComplexPolygonCell(index, box); |
490 | |
491 | auto first_half_planes = findHalfPlanes(box, first); |
492 | auto second_half_planes = findHalfPlanes(box, second); |
493 | |
494 | if (first_half_planes.empty()) |
495 | addCell(index, box, first); |
496 | else if (second_half_planes.empty()) |
497 | addCell(index, box, second); |
498 | else if (first_half_planes.size() == 1 && second_half_planes.size() == 1) |
499 | { |
500 | cells[index].type = CellType::pairOfLinesDifferentPolygons; |
501 | cells[index].half_planes[0] = first_half_planes[0]; |
502 | cells[index].half_planes[1] = second_half_planes[0]; |
503 | } |
504 | else |
505 | addComplexPolygonCell(index, box); |
506 | } |
507 | |
508 | |
509 | template <typename Strategy, typename CoordinateType = Float32> |
510 | class PointInPolygon |
511 | { |
512 | public: |
513 | using Point = boost::geometry::model::d2::point_xy<CoordinateType>; |
514 | /// Counter-Clockwise ordering. |
515 | using Polygon = boost::geometry::model::polygon<Point, false>; |
516 | using Box = boost::geometry::model::box<Point>; |
517 | |
518 | explicit PointInPolygon(const Polygon & polygon_) : polygon(polygon_) {} |
519 | |
520 | void init() |
521 | { |
522 | boost::geometry::envelope(polygon, box); |
523 | |
524 | const Point & min_corner = box.min_corner(); |
525 | const Point & max_corner = box.max_corner(); |
526 | |
527 | if (min_corner.x() == max_corner.x() || min_corner.y() == max_corner.y()) |
528 | has_empty_bound = true; |
529 | } |
530 | |
531 | bool hasEmptyBound() const { return has_empty_bound; } |
532 | |
533 | inline bool ALWAYS_INLINE contains(CoordinateType x, CoordinateType y) |
534 | { |
535 | Point point(x, y); |
536 | |
537 | if (!boost::geometry::within(point, box)) |
538 | return false; |
539 | |
540 | return boost::geometry::covered_by(point, polygon, strategy); |
541 | } |
542 | |
543 | UInt64 getAllocatedBytes() const { return sizeof(*this); } |
544 | |
545 | private: |
546 | const Polygon & polygon; |
547 | Box box; |
548 | bool has_empty_bound = false; |
549 | Strategy strategy; |
550 | }; |
551 | |
552 | |
553 | /// Algorithms. |
554 | |
555 | template <typename T, typename U, typename PointInPolygonImpl> |
556 | ColumnPtr pointInPolygon(const ColumnVector<T> & x, const ColumnVector<U> & y, PointInPolygonImpl && impl) |
557 | { |
558 | auto size = x.size(); |
559 | |
560 | impl.init(); |
561 | |
562 | if (impl.hasEmptyBound()) |
563 | { |
564 | return ColumnVector<UInt8>::create(size, 0); |
565 | } |
566 | |
567 | auto result = ColumnVector<UInt8>::create(size); |
568 | auto & data = result->getData(); |
569 | |
570 | const auto & x_data = x.getData(); |
571 | const auto & y_data = y.getData(); |
572 | |
573 | for (auto i : ext::range(0, size)) |
574 | { |
575 | data[i] = static_cast<UInt8>(impl.contains(x_data[i], y_data[i])); |
576 | } |
577 | |
578 | return result; |
579 | } |
580 | |
581 | template <typename ... Types> |
582 | struct CallPointInPolygon; |
583 | |
584 | template <typename Type, typename ... Types> |
585 | struct CallPointInPolygon<Type, Types ...> |
586 | { |
587 | template <typename T, typename PointInPolygonImpl> |
588 | static ColumnPtr call(const ColumnVector<T> & x, const IColumn & y, PointInPolygonImpl && impl) |
589 | { |
590 | if (auto column = typeid_cast<const ColumnVector<Type> *>(&y)) |
591 | return pointInPolygon(x, *column, impl); |
592 | return CallPointInPolygon<Types ...>::template call<T>(x, y, impl); |
593 | } |
594 | |
595 | template <typename PointInPolygonImpl> |
596 | static ColumnPtr call(const IColumn & x, const IColumn & y, PointInPolygonImpl && impl) |
597 | { |
598 | using Impl = typename ApplyTypeListForClass<::DB::GeoUtils::CallPointInPolygon, TypeListNativeNumbers>::Type; |
599 | if (auto column = typeid_cast<const ColumnVector<Type> *>(&x)) |
600 | return Impl::template call<Type>(*column, y, impl); |
601 | return CallPointInPolygon<Types ...>::call(x, y, impl); |
602 | } |
603 | }; |
604 | |
605 | template <> |
606 | struct CallPointInPolygon<> |
607 | { |
608 | template <typename T, typename PointInPolygonImpl> |
609 | static ColumnPtr call(const ColumnVector<T> &, const IColumn & y, PointInPolygonImpl &&) |
610 | { |
611 | throw Exception(std::string("Unknown numeric column type: " ) + demangle(typeid(y).name()), ErrorCodes::LOGICAL_ERROR); |
612 | } |
613 | |
614 | template <typename PointInPolygonImpl> |
615 | static ColumnPtr call(const IColumn & x, const IColumn &, PointInPolygonImpl &&) |
616 | { |
617 | throw Exception(std::string("Unknown numeric column type: " ) + demangle(typeid(x).name()), ErrorCodes::LOGICAL_ERROR); |
618 | } |
619 | }; |
620 | |
621 | template <typename PointInPolygonImpl> |
622 | ColumnPtr pointInPolygon(const IColumn & x, const IColumn & y, PointInPolygonImpl && impl) |
623 | { |
624 | using Impl = typename ApplyTypeListForClass<::DB::GeoUtils::CallPointInPolygon, TypeListNativeNumbers>::Type; |
625 | return Impl::call(x, y, impl); |
626 | } |
627 | |
628 | /// Total angle (signed) between neighbor vectors in linestring. Zero if linestring.size() < 2. |
629 | template <typename Linestring> |
630 | double calcLinestringRotation(const Linestring & points) |
631 | { |
632 | using Point = std::decay_t<decltype(*points.begin())>; |
633 | double rotation = 0; |
634 | |
635 | auto vecProduct = [](const Point & from, const Point & to) { return from.x() * to.y() - from.y() * to.x(); }; |
636 | auto scalarProduct = [](const Point & from, const Point & to) { return from.x() * to.x() + from.y() * to.y(); }; |
637 | auto getVector = [](const Point & from, const Point & to) -> Point |
638 | { |
639 | return Point(to.x() - from.x(), to.y() - from.y()); |
640 | }; |
641 | |
642 | for (auto it = points.begin(); it != points.end(); ++it) |
643 | { |
644 | if (it != points.begin()) |
645 | { |
646 | auto prev = std::prev(it); |
647 | auto next = std::next(it); |
648 | |
649 | if (next == points.end()) |
650 | next = std::next(points.begin()); |
651 | |
652 | Point from = getVector(*prev, *it); |
653 | Point to = getVector(*it, *next); |
654 | |
655 | auto vec_prod = vecProduct(from, to); |
656 | auto scalar_prod = scalarProduct(from, to); |
657 | auto ang = std::atan2(vec_prod, scalar_prod); |
658 | rotation += ang; |
659 | } |
660 | } |
661 | |
662 | return rotation; |
663 | } |
664 | |
665 | /// Make inner linestring counter-clockwise and outers clockwise oriented. |
666 | template <typename Polygon> |
667 | void normalizePolygon(Polygon && polygon) |
668 | { |
669 | auto & outer = polygon.outer(); |
670 | if (calcLinestringRotation(outer) < 0) |
671 | std::reverse(outer.begin(), outer.end()); |
672 | |
673 | auto & inners = polygon.inners(); |
674 | for (auto & inner : inners) |
675 | if (calcLinestringRotation(inner) > 0) |
676 | std::reverse(inner.begin(), inner.end()); |
677 | } |
678 | |
679 | |
680 | template <typename Polygon> |
681 | std::string serialize(Polygon && polygon) |
682 | { |
683 | std::string result; |
684 | |
685 | { |
686 | WriteBufferFromString buffer(result); |
687 | |
688 | using RingType = typename std::decay_t<Polygon>::ring_type; |
689 | |
690 | auto serializeFloat = [&buffer](float value) { buffer.write(reinterpret_cast<char *>(&value), sizeof(value)); }; |
691 | auto serializeSize = [&buffer](size_t size) { buffer.write(reinterpret_cast<char *>(&size), sizeof(size)); }; |
692 | |
693 | auto serializeRing = [& serializeFloat, & serializeSize](const RingType & ring) |
694 | { |
695 | serializeSize(ring.size()); |
696 | for (const auto & point : ring) |
697 | { |
698 | serializeFloat(point.x()); |
699 | serializeFloat(point.y()); |
700 | } |
701 | }; |
702 | |
703 | serializeRing(polygon.outer()); |
704 | |
705 | const auto & inners = polygon.inners(); |
706 | serializeSize(inners.size()); |
707 | for (auto & inner : inners) |
708 | serializeRing(inner); |
709 | } |
710 | |
711 | return result; |
712 | } |
713 | |
714 | size_t geohashEncode(Float64 longitude, Float64 latitude, UInt8 precision, char * out); |
715 | |
716 | void geohashDecode(const char * encoded_string, size_t encoded_len, Float64 * longitude, Float64 * latitude); |
717 | |
718 | std::vector<std::pair<Float64, Float64>> geohashCoverBox(Float64 longitude_min, Float64 latitude_min, Float64 longitude_max, Float64 latitude_max, UInt8 precision, UInt32 max_items = 0); |
719 | |
720 | struct GeohashesInBoxPreparedArgs |
721 | { |
722 | UInt64 items_count = 0; |
723 | UInt8 precision = 0; |
724 | |
725 | Float64 longitude_min = 0.0; |
726 | Float64 latitude_min = 0.0; |
727 | Float64 longitude_max = 0.0; |
728 | Float64 latitude_max = 0.0; |
729 | |
730 | Float64 longitude_step = 0.0; |
731 | Float64 latitude_step = 0.0; |
732 | }; |
733 | |
734 | GeohashesInBoxPreparedArgs geohashesInBoxPrepare(const Float64 longitude_min, |
735 | const Float64 latitude_min, |
736 | Float64 longitude_max, |
737 | Float64 latitude_max, |
738 | UInt8 precision); |
739 | |
740 | UInt64 geohashesInBox(const GeohashesInBoxPreparedArgs & estimation, char * out); |
741 | |
742 | } /// GeoUtils |
743 | |
744 | } /// DB |
745 | |