From 5508a4f89c20686a19f233ef0a04b796d8a2cbaa Mon Sep 17 00:00:00 2001 From: =?utf8?q?Ra=C3=BAl=20Mar=C3=ADn=20Rodr=C3=ADguez?= Date: Tue, 5 Feb 2019 15:27:53 +0000 Subject: [PATCH] Integrate wagyu to validate MVT polygons Closes https://github.com/postgis/postgis/pull/356 Closes #4311 git-svn-id: http://svn.osgeo.org/postgis/trunk@17231 b70326c6-7e19-0410-871a-916f4a2858ee --- .github/codecov.yml | 2 + .gitignore | 3 + .travis.yml | 1 + GNUmakefile.in | 2 +- NEWS | 1 + ci/travis/run_wagyu.sh | 18 + configure.ac | 48 +- deps/Makefile.in | 50 + deps/wagyu/LICENSE.README | 56 + deps/wagyu/LICENSE.geometry | 13 + deps/wagyu/LICENSE.wagyu | 37 + deps/wagyu/Makefile.in | 102 ++ deps/wagyu/README.md | 27 + deps/wagyu/include/mapbox/feature.hpp | 113 ++ deps/wagyu/include/mapbox/geometry.hpp | 12 + deps/wagyu/include/mapbox/geometry/box.hpp | 36 + deps/wagyu/include/mapbox/geometry/empty.hpp | 18 + .../include/mapbox/geometry/envelope.hpp | 33 + .../mapbox/geometry/for_each_point.hpp | 51 + .../include/mapbox/geometry/geometry.hpp | 56 + .../include/mapbox/geometry/line_string.hpp | 28 + .../mapbox/geometry/multi_line_string.hpp | 28 + .../include/mapbox/geometry/multi_point.hpp | 28 + .../include/mapbox/geometry/multi_polygon.hpp | 28 + deps/wagyu/include/mapbox/geometry/point.hpp | 42 + .../mapbox/geometry/point_arithmetic.hpp | 119 ++ .../wagyu/include/mapbox/geometry/polygon.hpp | 45 + .../geometry/wagyu/active_bound_list.hpp | 396 +++++ .../include/mapbox/geometry/wagyu/bound.hpp | 99 ++ .../mapbox/geometry/wagyu/bubble_sort.hpp | 28 + .../mapbox/geometry/wagyu/build_edges.hpp | 183 +++ .../wagyu/build_local_minima_list.hpp | 26 + .../mapbox/geometry/wagyu/build_result.hpp | 68 + .../include/mapbox/geometry/wagyu/config.hpp | 50 + .../include/mapbox/geometry/wagyu/edge.hpp | 120 ++ .../mapbox/geometry/wagyu/interrupt.hpp | 50 + .../mapbox/geometry/wagyu/intersect.hpp | 70 + .../mapbox/geometry/wagyu/intersect_util.hpp | 365 +++++ .../mapbox/geometry/wagyu/local_minimum.hpp | 117 ++ .../geometry/wagyu/local_minimum_util.hpp | 314 ++++ .../include/mapbox/geometry/wagyu/point.hpp | 110 ++ .../geometry/wagyu/process_horizontal.hpp | 256 ++++ .../mapbox/geometry/wagyu/process_maxima.hpp | 123 ++ .../mapbox/geometry/wagyu/quick_clip.hpp | 139 ++ .../include/mapbox/geometry/wagyu/ring.hpp | 633 ++++++++ .../mapbox/geometry/wagyu/ring_util.hpp | 832 ++++++++++ .../mapbox/geometry/wagyu/scanbeam.hpp | 45 + .../mapbox/geometry/wagyu/snap_rounding.hpp | 191 +++ .../geometry/wagyu/topology_correction.hpp | 1347 +++++++++++++++++ .../include/mapbox/geometry/wagyu/util.hpp | 89 ++ .../include/mapbox/geometry/wagyu/vatti.hpp | 64 + .../include/mapbox/geometry/wagyu/wagyu.hpp | 145 ++ deps/wagyu/include/mapbox/geometry_io.hpp | 98 ++ deps/wagyu/lwgeom_wagyu.cpp | 296 ++++ deps/wagyu/lwgeom_wagyu.h | 66 + doc/installation.xml | 10 + doc/reference_management.xml | 54 +- doc/reference_output.xml | 16 +- macros/ax_cxx_compile_stdcxx.m4 | 948 ++++++++++++ postgis/Makefile.in | 16 +- postgis/mvt.c | 48 +- postgis/postgis.sql.in | 10 + postgis/postgis_libprotobuf.c | 16 + postgis/postgis_module.c | 8 + postgis_config.h.in | 3 + regress/core/mvt.sql | 38 +- regress/core/mvt_expected | 4 +- utils/postgis_restore.pl.in | 1 + 68 files changed, 8464 insertions(+), 25 deletions(-) create mode 100644 ci/travis/run_wagyu.sh create mode 100644 deps/Makefile.in create mode 100644 deps/wagyu/LICENSE.README create mode 100644 deps/wagyu/LICENSE.geometry create mode 100644 deps/wagyu/LICENSE.wagyu create mode 100644 deps/wagyu/Makefile.in create mode 100644 deps/wagyu/README.md create mode 100644 deps/wagyu/include/mapbox/feature.hpp create mode 100644 deps/wagyu/include/mapbox/geometry.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/box.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/empty.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/envelope.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/for_each_point.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/geometry.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/line_string.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/multi_line_string.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/multi_point.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/multi_polygon.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/point.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/point_arithmetic.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/polygon.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/active_bound_list.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/bound.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/bubble_sort.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/build_edges.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/build_local_minima_list.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/build_result.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/config.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/edge.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/interrupt.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/intersect.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/intersect_util.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/local_minimum.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/local_minimum_util.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/point.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/process_horizontal.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/process_maxima.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/quick_clip.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/ring.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/ring_util.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/scanbeam.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/snap_rounding.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/topology_correction.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/util.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/vatti.hpp create mode 100644 deps/wagyu/include/mapbox/geometry/wagyu/wagyu.hpp create mode 100644 deps/wagyu/include/mapbox/geometry_io.hpp create mode 100644 deps/wagyu/lwgeom_wagyu.cpp create mode 100644 deps/wagyu/lwgeom_wagyu.h create mode 100644 macros/ax_cxx_compile_stdcxx.m4 diff --git a/.github/codecov.yml b/.github/codecov.yml index f440b76e7..2e2529ecb 100644 --- a/.github/codecov.yml +++ b/.github/codecov.yml @@ -7,3 +7,5 @@ coverage: default: threshold: 0.2% comment: false +ignore: + - "deps/wagyu/include/**/*" # Ignore wagyu library files diff --git a/.gitignore b/.gitignore index 4905e6305..763d85069 100644 --- a/.gitignore +++ b/.gitignore @@ -182,3 +182,6 @@ postgis/uninstall_sfcgal.sql # LLVM JIT *.bc *.ll + +deps/Makefile +deps/wagyu/Makefile diff --git a/.travis.yml b/.travis.yml index 4c402b6da..698869b8f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -10,6 +10,7 @@ env: - tag=pg11-geos37-gdal24-proj52 mode=coverage - tag=pg11-geos37-gdal24-proj52 mode=usan_gcc - tag=pg11-geos37-gdal24-proj52 mode=usan_clang + - tag=pg11-geos37-gdal24-proj52 mode=wagyu - tag=pg11-geos37-gdal24-proj52 mode=tests - tag=pg10-geos36-gdal23-proj49 mode=tests - tag=pg96-geos36-gdal22-proj49 mode=tests diff --git a/GNUmakefile.in b/GNUmakefile.in index 3469a31af..f05782be1 100644 --- a/GNUmakefile.in +++ b/GNUmakefile.in @@ -6,7 +6,7 @@ SUBDIRS = liblwgeom ifeq (@LIBLWGEOM_ONLY@,no) -SUBDIRS += libpgcommon postgis regress @RASTER@ @TOPOLOGY@ loader utils doc @EXTENSIONS@ +SUBDIRS += libpgcommon postgis regress @RASTER@ @TOPOLOGY@ loader utils doc @EXTENSIONS@ @DEPS_SUBDIR@ endif POSTGIS_MAJOR_VERSION=@POSTGIS_MAJOR_VERSION@ diff --git a/NEWS b/NEWS index 3b233d137..835b3632e 100644 --- a/NEWS +++ b/NEWS @@ -26,6 +26,7 @@ PostGIS 3.0.0 within, equals (Esteban Zimányi and Arthur Lesuisse from Université Libre de Bruxelles (ULB), Darafei Praliaskouski) - #4171, ST_3DLineInterpolatePoint (Julien Cabieces, Vincent Mora) + - #4311, Introduce `--with-wagyu` as an option for MVT polygons (Raúl Marín) * Enhancements and fixes * - #4153, ST_Segmentize now splits segments proportionally (Darafei diff --git a/ci/travis/run_wagyu.sh b/ci/travis/run_wagyu.sh new file mode 100644 index 000000000..e76fbd03f --- /dev/null +++ b/ci/travis/run_wagyu.sh @@ -0,0 +1,18 @@ +#!/usr/bin/env bash +set -e + +WARNINGS="-Werror -Wall -Wextra -Wformat -Werror=format-security" +WARNINGS_DISABLED="-Wno-unused-parameter -Wno-implicit-fallthrough -Wno-cast-function-type" + +# Build with coverage +CFLAGS="-g -O0 --coverage -mtune=generic -fno-omit-frame-pointer ${WARNINGS} ${WARNINGS_DISABLED}" +LDFLAGS="--coverage" + +/usr/local/pgsql/bin/pg_ctl -c -l /tmp/logfile start +./autogen.sh + +# Standard build +./configure --with-wagyu CFLAGS="${CFLAGS}" CXXFLAGS="${CFLAGS}" LDFLAGS="${LDFLAGS}" +bash ./ci/travis/logbt -- make -j check RUNTESTFLAGS=--verbose +curl -S -f https://codecov.io/bash -o .github/codecov.bash +bash .github/codecov.bash diff --git a/configure.ac b/configure.ac index 2583c4ce9..37f97fb0e 100644 --- a/configure.ac +++ b/configure.ac @@ -1421,6 +1421,46 @@ else fi +dnl =========================================================================== +dnl Deps folder +dnl =========================================================================== + +DEPS_MAKEFILE_LIST=" + deps/Makefile" + + +dnl =========================================================================== +dnl Wagyu +dnl =========================================================================== + +dnl Wagyu will only be necessary if protobuf is present to build MVTs + +HAVE_WAGYU=no +if test "x$HAVE_PROTOBUF" = "xyes"; then + AC_ARG_WITH( + [wagyu], + AC_HELP_STRING([--with-wagyu], [Use the wagyu library]), + [HAVE_WAGYU=yes], + [HAVE_WAGYU=no]) + + if test "x$HAVE_WAGYU" = "xyes"; then + AC_MSG_RESULT([WAGYU: Wagyu usage requested]) + + DEPS_SUBDIR="deps" + AC_SUBST([DEPS_SUBDIR]) + + WAGYU_LIB=libwagyu.la + AC_SUBST([WAGYU_LIB]) + + AC_PROG_CXX + AX_CXX_COMPILE_STDCXX(11, noext, mandatory) + + AC_DEFINE([HAVE_WAGYU], [1], [Define to 1 if wagyu is being built]) + AC_SUBST([HAVE_WAGYU]) + DEPS_MAKEFILE_LIST="$DEPS_MAKEFILE_LIST + deps/wagyu/Makefile" + fi +fi @@ -1471,7 +1511,8 @@ AC_OUTPUT([GNUmakefile doc/Makefile.comments doc/html/image_src/Makefile utils/Makefile - $RT_MAKEFILE_LIST]) + $RT_MAKEFILE_LIST + $DEPS_MAKEFILE_LIST]) dnl =========================================================================== dnl Display the configuration status information @@ -1482,6 +1523,9 @@ AC_MSG_RESULT([ PostGIS is now configured for ${host}]) AC_MSG_RESULT() AC_MSG_RESULT([ -------------- Compiler Info ------------- ]) AC_MSG_RESULT([ C compiler: ${CC} ${CFLAGS}]) +if test "x$HAVE_WAGYU" = "xyes"; then + AC_MSG_RESULT([ C++ compiler: ${CXX} ${CXXFLAGS}]) +fi AC_MSG_RESULT([ SQL preprocessor: ${SQLPP}]) AC_MSG_RESULT() AC_MSG_RESULT([ -------------- Additional Info ------------- ]) @@ -1520,6 +1564,8 @@ AC_MSG_RESULT([ JSON-C support: ${HAVE_JSON}]) AC_MSG_RESULT([ protobuf-c support: ${HAVE_PROTOBUF}]) AC_MSG_RESULT([ PCRE support: ${HAVE_PCRE}]) AC_MSG_RESULT([ Perl: ${PERL}]) +AC_MSG_RESULT([ Wagyu: ${HAVE_WAGYU}]) + AC_MSG_RESULT() AC_MSG_RESULT([ --------------- Extensions --------------- ]) if test "x$RASTER" = "xraster"; then diff --git a/deps/Makefile.in b/deps/Makefile.in new file mode 100644 index 000000000..25bdddf35 --- /dev/null +++ b/deps/Makefile.in @@ -0,0 +1,50 @@ +#/********************************************************************** +# * +# * PostGIS - Spatial Types for PostgreSQL +# * http://postgis.net +# * +# * PostGIS is free software: you can redistribute it and/or modify +# * it under the terms of the GNU General Public License as published by +# * the Free Software Foundation, either version 2 of the License, or +# * (at your option) any later version. +# * +# * PostGIS is distributed in the hope that it will be useful, +# * but WITHOUT ANY WARRANTY; without even the implied warranty of +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# * GNU General Public License for more details. +# * +# * You should have received a copy of the GNU General Public License +# * along with PostGIS. If not, see . +# * +# ********************************************************************** +# * +# * Copyright 2019 Raúl Marín +# * +# **********************************************************************/ + +CC=@CC@ +CXX=@CXX@ +CFLAGS=-I../liblwgeom @WARNFLAGS@ @CFLAGS@ @PICFLAGS@ +CXXFLAGS=-I../liblwgeom @WARNFLAGS@ @CXXFLAGS@ @PICFLAGS@ +LDFLAGS = @LDFLAGS@ +top_builddir = @top_builddir@ +libdir = @libdir@ +LIBTOOL = @LIBTOOL@ + +all: @WAGYU_LIB@ + +@WAGYU_LIB@: + $(MAKE) -C wagyu $@ + +install: + +uninstall: + +check: + +clean: + $(MAKE) -C wagyu $@ + +distclean: clean + $(MAKE) -C wagyu $@ + rm -f Makefile diff --git a/deps/wagyu/LICENSE.README b/deps/wagyu/LICENSE.README new file mode 100644 index 000000000..23ce942b6 --- /dev/null +++ b/deps/wagyu/LICENSE.README @@ -0,0 +1,56 @@ +# LICENSE.geometry applies to: + +include/ +└── mapbox + ├── feature.hpp + ├── geometry + │   ├── box.hpp + │   ├── empty.hpp + │   ├── envelope.hpp + │   ├── for_each_point.hpp + │   ├── geometry.hpp + │   ├── line_string.hpp + │   ├── multi_line_string.hpp + │   ├── multi_point.hpp + │   ├── multi_polygon.hpp + │   ├── point_arithmetic.hpp + │   ├── point.hpp + │   └── polygon.hpp + ├── geometry.hpp + └── geometry_io.hpp + +2 directories, 15 files + +# LICENSE.wagyu applies to: + +include/ +└── mapbox + └── geometry + └── wagyu + ├── active_bound_list.hpp + ├── bound.hpp + ├── bubble_sort.hpp + ├── build_edges.hpp + ├── build_local_minima_list.hpp + ├── build_result.hpp + ├── config.hpp + ├── edge.hpp + ├── interrupt.hpp + ├── intersect.hpp + ├── intersect_util.hpp + ├── local_minimum.hpp + ├── local_minimum_util.hpp + ├── point.hpp + ├── process_horizontal.hpp + ├── process_maxima.hpp + ├── quick_clip.hpp + ├── ring.hpp + ├── ring_util.hpp + ├── scanbeam.hpp + ├── snap_rounding.hpp + ├── topology_correction.hpp + ├── util.hpp + ├── vatti.hpp + └── wagyu.hpp + +3 directories, 24 files diff --git a/deps/wagyu/LICENSE.geometry b/deps/wagyu/LICENSE.geometry new file mode 100644 index 000000000..8f6a86bd0 --- /dev/null +++ b/deps/wagyu/LICENSE.geometry @@ -0,0 +1,13 @@ +Copyright (c) 2016, Mapbox + +Permission to use, copy, modify, and/or distribute this software for any +purpose with or without fee is hereby granted, provided that the above +copyright notice and this permission notice appear in all copies. + +THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES +WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR +ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES +WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN +ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF +OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. \ No newline at end of file diff --git a/deps/wagyu/LICENSE.wagyu b/deps/wagyu/LICENSE.wagyu new file mode 100644 index 000000000..f3eb477c3 --- /dev/null +++ b/deps/wagyu/LICENSE.wagyu @@ -0,0 +1,37 @@ +Parts of the code in the Wagyu Library are derived from the version of the +Clipper Library by Angus Johnson listed below. + +Author : Angus Johnson +Version : 6.4.0 +Date : 2 July 2015 +Website : http://www.angusj.com + +Copyright for portions of the derived code in the Wagyu library are held +by Angus Johnson, 2010-2015. All other copyright for the Wagyu Library are held by +Mapbox, 2016. This code is published in accordance with, and retains the same license +as the Clipper Library by Angus Johnson. + +Copyright (c) 2010-2015, Angus Johnson +Copyright (c) 2016, Mapbox + +Permission is hereby granted, free of charge, to any person or organization +obtaining a copy of the software and accompanying documentation covered by +this license (the "Software") to use, reproduce, display, distribute, +execute, and transmit the Software, and to prepare derivative works of the +Software, and to permit third-parties to whom the Software is furnished to +do so, all subject to the following: + +The copyright notices in the Software and this entire statement, including +the above license grant, this restriction and the following disclaimer, +must be included in all copies of the Software, in whole or in part, and +all derivative works of the Software, unless such copies or derivative +works are solely in the form of machine-executable object code generated by +a source language processor. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT +SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE +FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, +ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. diff --git a/deps/wagyu/Makefile.in b/deps/wagyu/Makefile.in new file mode 100644 index 000000000..1e3d1a4bc --- /dev/null +++ b/deps/wagyu/Makefile.in @@ -0,0 +1,102 @@ +#/********************************************************************** +# * +# * PostGIS - Spatial Types for PostgreSQL +# * http://postgis.net +# * +# * PostGIS is free software: you can redistribute it and/or modify +# * it under the terms of the GNU General Public License as published by +# * the Free Software Foundation, either version 2 of the License, or +# * (at your option) any later version. +# * +# * PostGIS is distributed in the hope that it will be useful, +# * but WITHOUT ANY WARRANTY; without even the implied warranty of +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# * GNU General Public License for more details. +# * +# * You should have received a copy of the GNU General Public License +# * along with PostGIS. If not, see . +# * +# ********************************************************************** +# * +# * Copyright 2019 Raúl Marín +# * +# **********************************************************************/ + +CXX = @CXX@ +CXXFLAGS =-I../../liblwgeom -Iinclude @CPPFLAGS@ @WARNFLAGS@ @CXXFLAGS@ @PICFLAGS@ +LDFLAGS = @LDFLAGS@ +top_builddir = @top_builddir@ +libdir = @libdir@ +LIBTOOL = @LIBTOOL@ + +WAGYU_OBJS = \ + lwgeom_wagyu.o + +WAGYU_HEADERS = \ + lwgeom_wagyu.h \ + include/mapbox/geometry/polygon.hpp \ + include/mapbox/geometry/point_arithmetic.hpp \ + include/mapbox/geometry/multi_polygon.hpp \ + include/mapbox/geometry/box.hpp \ + include/mapbox/geometry/for_each_point.hpp \ + include/mapbox/geometry/line_string.hpp \ + include/mapbox/geometry/wagyu/process_horizontal.hpp \ + include/mapbox/geometry/wagyu/wagyu.hpp \ + include/mapbox/geometry/wagyu/ring_util.hpp \ + include/mapbox/geometry/wagyu/bubble_sort.hpp \ + include/mapbox/geometry/wagyu/snap_rounding.hpp \ + include/mapbox/geometry/wagyu/bound.hpp \ + include/mapbox/geometry/wagyu/edge.hpp \ + include/mapbox/geometry/wagyu/build_result.hpp \ + include/mapbox/geometry/wagyu/scanbeam.hpp \ + include/mapbox/geometry/wagyu/ring.hpp \ + include/mapbox/geometry/wagyu/config.hpp \ + include/mapbox/geometry/wagyu/quick_clip.hpp \ + include/mapbox/geometry/wagyu/local_minimum_util.hpp \ + include/mapbox/geometry/wagyu/process_maxima.hpp \ + include/mapbox/geometry/wagyu/intersect.hpp \ + include/mapbox/geometry/wagyu/point.hpp \ + include/mapbox/geometry/wagyu/util.hpp \ + include/mapbox/geometry/wagyu/intersect_util.hpp \ + include/mapbox/geometry/wagyu/vatti.hpp \ + include/mapbox/geometry/wagyu/active_bound_list.hpp \ + include/mapbox/geometry/wagyu/build_local_minima_list.hpp \ + include/mapbox/geometry/wagyu/build_edges.hpp \ + include/mapbox/geometry/wagyu/local_minimum.hpp \ + include/mapbox/geometry/wagyu/topology_correction.hpp \ + include/mapbox/geometry/wagyu/interrupt.hpp \ + include/mapbox/geometry/geometry.hpp \ + include/mapbox/geometry/point.hpp \ + include/mapbox/geometry/empty.hpp \ + include/mapbox/geometry/multi_line_string.hpp \ + include/mapbox/geometry/envelope.hpp \ + include/mapbox/geometry/multi_point.hpp \ + include/mapbox/feature.hpp \ + include/mapbox/geometry_io.hpp \ + include/mapbox/geometry.hpp + +all: @WAGYU_LIB@ + +@WAGYU_LIB@: $(WAGYU_OBJS) + ar rs @WAGYU_LIB@ $(WAGYU_OBJS) + +$(WAGYU_OBJS): %.o: %.cpp ../../liblwgeom/liblwgeom.h $(WAGYU_HEADERS) + $(CXX) $(CXXFLAGS) -c -o $@ $< + +../../liblwgeom/liblwgeom.h: + $(MAKE) -C ../../liblwgeom liblwgeom.h + +clean: + rm -f @WAGYU_LIB@ $(WAGYU_OBJS) + +distclean: clean + rm -f Makefile + + +install: + +uninstall: + +check: + +.PHONY: clean distclean install uninstall check diff --git a/deps/wagyu/README.md b/deps/wagyu/README.md new file mode 100644 index 000000000..94a0edad1 --- /dev/null +++ b/deps/wagyu/README.md @@ -0,0 +1,27 @@ +# Wagyu + +Wagyu is a library that performs polygon clipping using the [Vatti Clipping Algorithm](https://en.wikipedia.org/wiki/Vatti_clipping_algorithm). It is based on the [Clipper Library](http://www.angusj.com/delphi/clipper.php) but it adds an extra phase to correct the topology and guarantee validation accoding to the OGC specification. + +Wagyu is a header only library but depends on [Mapbox Geometry](https://github.com/mapbox/geometry.hpp). It requires a C++ compiler that supports at least C++11. + +# liblwgeom bindings + +To be able to use Wagyu in Postgis, several functions have been added as a bridge between liblwgeom's C code, and the C++ header only library: + - `lwgeom_wagyu_clip_by_box`: Main function to clip and force validation over a polygon. + - `libwagyu_version`: Returns a static string with the version of the wagyu library. + - `lwgeom_wagyu_interruptRequest`: Function to request the interruption of the wagyu library. + +It is integrated in the project as an static library inside postgis.so, so it doesn't require an extra dependency at runtime besides `libstdc++` and `libm` which were a requisite. + +# Main considerations about the library + +- It works only with POLYGONTYPE or MULTIPOLYGONTYPE type geometries. Anything else will be dropped. +- It's currently setup to use `int32_t` values internally for extra speed. It could be changed to use to `int64_t` or to `double` to match liblwgeom but, as it's only used by MVT, this isn't necessary. +- The library is clipping the geometry to the bounding box passed. It also supports `Intersection` of 2 geometries, `Union`, `Difference` or `XOR` of 2 geometries but those functionalities aren't exposed. +- The library outputs the geometry in the correct winding order for **MVT**, which is the opposite of OGC as the vertical order is switched. + + +# Dependency changelog + + - 2019-02-05 - [Wagyu] Library extraction from https://github.com/mapbox/wagyu + - 2019-02-05 - [geometry.hpp] Library extraction from https://github.com/mapbox/geometry.hpp diff --git a/deps/wagyu/include/mapbox/feature.hpp b/deps/wagyu/include/mapbox/feature.hpp new file mode 100644 index 000000000..8d03f4674 --- /dev/null +++ b/deps/wagyu/include/mapbox/feature.hpp @@ -0,0 +1,113 @@ +#pragma once + +#include + +#include + +#include +#include +#include +#include + +namespace mapbox { +namespace feature { + +struct value; + +struct null_value_t +{ +}; + +constexpr bool operator==(const null_value_t&, const null_value_t&) { return true; } +constexpr bool operator!=(const null_value_t&, const null_value_t&) { return false; } +constexpr bool operator<(const null_value_t&, const null_value_t&) { return false; } + +constexpr null_value_t null_value = null_value_t(); + +// Multiple numeric types (uint64_t, int64_t, double) are present in order to support +// the widest possible range of JSON numbers, which do not have a maximum range. +// Implementations that produce `value`s should use that order for type preference, +// using uint64_t for positive integers, int64_t for negative integers, and double +// for non-integers and integers outside the range of 64 bits. +using value_base = mapbox::util::variant>, + mapbox::util::recursive_wrapper>>; + +struct value : value_base +{ + using value_base::value_base; +}; + +using property_map = std::unordered_map; + +// The same considerations and requirement for numeric types apply as for `value_base`. +using identifier = mapbox::util::variant; + +template +struct feature +{ + using coordinate_type = T; + using geometry_type = mapbox::geometry::geometry; // Fully qualified to avoid GCC -fpermissive error. + + geometry_type geometry; + property_map properties; + identifier id; + + feature() + : geometry(), + properties(), + id() {} + feature(geometry_type const& geom_) + : geometry(geom_), + properties(), + id() {} + feature(geometry_type&& geom_) + : geometry(std::move(geom_)), + properties(), + id() {} + feature(geometry_type const& geom_, property_map const& prop_) + : geometry(geom_), properties(prop_), id() {} + feature(geometry_type&& geom_, property_map&& prop_) + : geometry(std::move(geom_)), + properties(std::move(prop_)), + id() {} + feature(geometry_type const& geom_, property_map const& prop_, identifier const& id_) + : geometry(geom_), + properties(prop_), + id(id_) {} + feature(geometry_type&& geom_, property_map&& prop_, identifier&& id_) + : geometry(std::move(geom_)), + properties(std::move(prop_)), + id(std::move(id_)) {} +}; + +template +constexpr bool operator==(feature const& lhs, feature const& rhs) +{ + return lhs.id == rhs.id && lhs.geometry == rhs.geometry && lhs.properties == rhs.properties; +} + +template +constexpr bool operator!=(feature const& lhs, feature const& rhs) +{ + return !(lhs == rhs); +} + +template class Cont = std::vector> +struct feature_collection : Cont> +{ + using coordinate_type = T; + using feature_type = feature; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + feature_collection(Args&&... args) : container_type(std::forward(args)...) + { + } + feature_collection(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +} // namespace feature +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry.hpp b/deps/wagyu/include/mapbox/geometry.hpp new file mode 100644 index 000000000..9caecd9a1 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include diff --git a/deps/wagyu/include/mapbox/geometry/box.hpp b/deps/wagyu/include/mapbox/geometry/box.hpp new file mode 100644 index 000000000..bb3ea5701 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/box.hpp @@ -0,0 +1,36 @@ +#pragma once + +#include + +namespace mapbox { +namespace geometry { + +template +struct box +{ + using coordinate_type = T; + using point_type = point; + + constexpr box(point_type const& min_, point_type const& max_) + : min(min_), max(max_) + { + } + + point_type min; + point_type max; +}; + +template +constexpr bool operator==(box const& lhs, box const& rhs) +{ + return lhs.min == rhs.min && lhs.max == rhs.max; +} + +template +constexpr bool operator!=(box const& lhs, box const& rhs) +{ + return lhs.min != rhs.min || lhs.max != rhs.max; +} + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/empty.hpp b/deps/wagyu/include/mapbox/geometry/empty.hpp new file mode 100644 index 000000000..0ea446d8b --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/empty.hpp @@ -0,0 +1,18 @@ +#pragma once + +namespace mapbox { +namespace geometry { + +struct empty +{ +}; // this Geometry type represents the empty point set, ∅, for the coordinate space (OGC Simple Features). + +constexpr bool operator==(empty, empty) { return true; } +constexpr bool operator!=(empty, empty) { return false; } +constexpr bool operator<(empty, empty) { return false; } +constexpr bool operator>(empty, empty) { return false; } +constexpr bool operator<=(empty, empty) { return true; } +constexpr bool operator>=(empty, empty) { return true; } + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/envelope.hpp b/deps/wagyu/include/mapbox/geometry/envelope.hpp new file mode 100644 index 000000000..e1f722fc5 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/envelope.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include +#include + +#include + +namespace mapbox { +namespace geometry { + +template +box envelope(G const& geometry) +{ + using limits = std::numeric_limits; + + T min_t = limits::has_infinity ? -limits::infinity() : limits::min(); + T max_t = limits::has_infinity ? limits::infinity() : limits::max(); + + point min(max_t, max_t); + point max(min_t, min_t); + + for_each_point(geometry, [&](point const& point) { + if (min.x > point.x) min.x = point.x; + if (min.y > point.y) min.y = point.y; + if (max.x < point.x) max.x = point.x; + if (max.y < point.y) max.y = point.y; + }); + + return box(min, max); +} + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/for_each_point.hpp b/deps/wagyu/include/mapbox/geometry/for_each_point.hpp new file mode 100644 index 000000000..d31b484ff --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/for_each_point.hpp @@ -0,0 +1,51 @@ +#pragma once + +#include + +namespace mapbox { +namespace geometry { + +template +void for_each_point(mapbox::geometry::empty const&, F&&) +{ +} + +template +auto for_each_point(Point&& point, F&& f) + -> decltype(point.x, point.y, void()) +{ + f(std::forward(point)); +} + +template +auto for_each_point(Container&& container, F&& f) + -> decltype(container.begin(), container.end(), void()); + +template +void for_each_point(mapbox::util::variant const& geom, F&& f) +{ + mapbox::util::variant::visit(geom, [&](auto const& g) { + for_each_point(g, f); + }); +} + +template +void for_each_point(mapbox::util::variant& geom, F&& f) +{ + mapbox::util::variant::visit(geom, [&](auto& g) { + for_each_point(g, f); + }); +} + +template +auto for_each_point(Container&& container, F&& f) + -> decltype(container.begin(), container.end(), void()) +{ + for (auto& e : container) + { + for_each_point(e, f); + } +} + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/geometry.hpp b/deps/wagyu/include/mapbox/geometry/geometry.hpp new file mode 100644 index 000000000..fa1e7c31c --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/geometry.hpp @@ -0,0 +1,56 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include + +// stl +#include + +namespace mapbox { +namespace geometry { + +template class Cont = std::vector> +struct geometry_collection; + +template class Cont = std::vector> +using geometry_base = mapbox::util::variant, + line_string, + polygon, + multi_point, + multi_line_string, + multi_polygon, + geometry_collection>; + +template class Cont = std::vector> +struct geometry : geometry_base +{ + using coordinate_type = T; + using geometry_base::geometry_base; +}; + +template class Cont> +struct geometry_collection : Cont> +{ + using coordinate_type = T; + using geometry_type = geometry; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + geometry_collection(Args&&... args) : container_type(std::forward(args)...) + { + } + geometry_collection(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/line_string.hpp b/deps/wagyu/include/mapbox/geometry/line_string.hpp new file mode 100644 index 000000000..a015f71d1 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/line_string.hpp @@ -0,0 +1,28 @@ +#pragma once + +// mapbox +#include +// stl +#include + +namespace mapbox { +namespace geometry { + +template class Cont = std::vector> +struct line_string : Cont> +{ + using coordinate_type = T; + using point_type = point; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + line_string(Args&&... args) : container_type(std::forward(args)...) + { + } + line_string(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/multi_line_string.hpp b/deps/wagyu/include/mapbox/geometry/multi_line_string.hpp new file mode 100644 index 000000000..b528fa61b --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/multi_line_string.hpp @@ -0,0 +1,28 @@ +#pragma once + +// mapbox +#include +// stl +#include + +namespace mapbox { +namespace geometry { + +template class Cont = std::vector> +struct multi_line_string : Cont> +{ + using coordinate_type = T; + using line_string_type = line_string; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + multi_line_string(Args&&... args) : container_type(std::forward(args)...) + { + } + multi_line_string(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/multi_point.hpp b/deps/wagyu/include/mapbox/geometry/multi_point.hpp new file mode 100644 index 000000000..060927bf2 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/multi_point.hpp @@ -0,0 +1,28 @@ +#pragma once + +// mapbox +#include +// stl +#include + +namespace mapbox { +namespace geometry { + +template class Cont = std::vector> +struct multi_point : Cont> +{ + using coordinate_type = T; + using point_type = point; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + multi_point(Args&&... args) : container_type(std::forward(args)...) + { + } + multi_point(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/multi_polygon.hpp b/deps/wagyu/include/mapbox/geometry/multi_polygon.hpp new file mode 100644 index 000000000..9546860ed --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/multi_polygon.hpp @@ -0,0 +1,28 @@ +#pragma once + +// mapbox +#include +// stl +#include + +namespace mapbox { +namespace geometry { + +template class Cont = std::vector> +struct multi_polygon : Cont> +{ + using coordinate_type = T; + using polygon_type = polygon; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + multi_polygon(Args&&... args) : container_type(std::forward(args)...) + { + } + multi_polygon(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/point.hpp b/deps/wagyu/include/mapbox/geometry/point.hpp new file mode 100644 index 000000000..da8d67733 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/point.hpp @@ -0,0 +1,42 @@ +#pragma once + +namespace mapbox { +namespace geometry { + +template +struct point +{ + using coordinate_type = T; + + constexpr point() + : x(), y() + { + } + constexpr point(T x_, T y_) + : x(x_), y(y_) + { + } + + T x; + T y; +}; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wfloat-equal" + +template +constexpr bool operator==(point const& lhs, point const& rhs) +{ + return lhs.x == rhs.x && lhs.y == rhs.y; +} + +#pragma GCC diagnostic pop + +template +constexpr bool operator!=(point const& lhs, point const& rhs) +{ + return !(lhs == rhs); +} + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/point_arithmetic.hpp b/deps/wagyu/include/mapbox/geometry/point_arithmetic.hpp new file mode 100644 index 000000000..0c4c63278 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/point_arithmetic.hpp @@ -0,0 +1,119 @@ +#pragma once + +namespace mapbox { +namespace geometry { + +template +point operator+(point const& lhs, point const& rhs) +{ + return point(lhs.x + rhs.x, lhs.y + rhs.y); +} + +template +point operator+(point const& lhs, T const& rhs) +{ + return point(lhs.x + rhs, lhs.y + rhs); +} + +template +point operator-(point const& lhs, point const& rhs) +{ + return point(lhs.x - rhs.x, lhs.y - rhs.y); +} + +template +point operator-(point const& lhs, T const& rhs) +{ + return point(lhs.x - rhs, lhs.y - rhs); +} + +template +point operator*(point const& lhs, point const& rhs) +{ + return point(lhs.x * rhs.x, lhs.y * rhs.y); +} + +template +point operator*(point const& lhs, T const& rhs) +{ + return point(lhs.x * rhs, lhs.y * rhs); +} + +template +point operator/(point const& lhs, point const& rhs) +{ + return point(lhs.x / rhs.x, lhs.y / rhs.y); +} + +template +point operator/(point const& lhs, T const& rhs) +{ + return point(lhs.x / rhs, lhs.y / rhs); +} + +template +point& operator+=(point& lhs, point const& rhs) +{ + lhs.x += rhs.x; + lhs.y += rhs.y; + return lhs; +} + +template +point& operator+=(point& lhs, T const& rhs) +{ + lhs.x += rhs; + lhs.y += rhs; + return lhs; +} + +template +point& operator-=(point& lhs, point const& rhs) +{ + lhs.x -= rhs.x; + lhs.y -= rhs.y; + return lhs; +} + +template +point& operator-=(point& lhs, T const& rhs) +{ + lhs.x -= rhs; + lhs.y -= rhs; + return lhs; +} + +template +point& operator*=(point& lhs, point const& rhs) +{ + lhs.x *= rhs.x; + lhs.y *= rhs.y; + return lhs; +} + +template +point& operator*=(point& lhs, T const& rhs) +{ + lhs.x *= rhs; + lhs.y *= rhs; + return lhs; +} + +template +point& operator/=(point& lhs, point const& rhs) +{ + lhs.x /= rhs.x; + lhs.y /= rhs.y; + return lhs; +} + +template +point& operator/=(point& lhs, T const& rhs) +{ + lhs.x /= rhs; + lhs.y /= rhs; + return lhs; +} + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/polygon.hpp b/deps/wagyu/include/mapbox/geometry/polygon.hpp new file mode 100644 index 000000000..5cad7bbf5 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/polygon.hpp @@ -0,0 +1,45 @@ +#pragma once + +// mapbox +#include + +// stl +#include + +namespace mapbox { +namespace geometry { + +template class Cont = std::vector> +struct linear_ring : Cont> +{ + using coordinate_type = T; + using point_type = point; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + linear_ring(Args&&... args) : container_type(std::forward(args)...) + { + } + linear_ring(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +template class Cont = std::vector> +struct polygon : Cont> +{ + using coordinate_type = T; + using linear_ring_type = linear_ring; + using container_type = Cont; + using size_type = typename container_type::size_type; + + template + polygon(Args&&... args) : container_type(std::forward(args)...) + { + } + polygon(std::initializer_list args) + : container_type(std::move(args)) {} +}; + +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/active_bound_list.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/active_bound_list.hpp new file mode 100644 index 000000000..dc113eea8 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/active_bound_list.hpp @@ -0,0 +1,396 @@ +#pragma once + +#ifdef DEBUG +#include +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +using active_bound_list = std::vector>; + +template +using active_bound_list_itr = typename active_bound_list::iterator; + +template +using active_bound_list_rev_itr = typename active_bound_list::reverse_iterator; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const active_bound_list& bnds) { + std::size_t c = 0; + for (auto const& bnd : bnds) { + out << "Index: " << c++ << std::endl; + out << *bnd; + } + return out; +} + +template +std::string output_edges(active_bound_list const& bnds) { + std::ostringstream out; + out << "["; + bool first = true; + for (auto const& bnd : bnds) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << bnd->current_edge->bot.x << "," << bnd->current_edge->bot.y << "],["; + out << bnd->current_edge->top.x << "," << bnd->current_edge->top.y << "]]"; + } + out << "]"; + return out.str(); +} + +#endif + +template +bool is_even_odd_fill_type(bound const& bound, fill_type subject_fill_type, fill_type clip_fill_type) { + if (bound.poly_type == polygon_type_subject) { + return subject_fill_type == fill_type_even_odd; + } else { + return clip_fill_type == fill_type_even_odd; + } +} + +template +bool is_even_odd_alt_fill_type(bound const& bound, fill_type subject_fill_type, fill_type clip_fill_type) { + if (bound.poly_type == polygon_type_subject) { + return clip_fill_type == fill_type_even_odd; + } else { + return subject_fill_type == fill_type_even_odd; + } +} + +template +struct bound_insert_location { + + bound const& bound2; + + bound_insert_location(bound const& b) : bound2(b) { + } + + bool operator()(bound_ptr const& b) { + auto const& bound1 = *b; + if (values_are_equal(bound2.current_x, bound1.current_x)) { + if (bound2.current_edge->top.y > bound1.current_edge->top.y) { + return static_cast(bound2.current_edge->top.x) < + get_current_x(*(bound1.current_edge), bound2.current_edge->top.y); + } else { + return static_cast(bound1.current_edge->top.x) > + get_current_x(*(bound2.current_edge), bound1.current_edge->top.y); + } + } else { + return bound2.current_x < bound1.current_x; + } + } +}; + +template +active_bound_list_itr insert_bound_into_ABL(bound& left, bound& right, active_bound_list& active_bounds) { + + auto itr = std::find_if(active_bounds.begin(), active_bounds.end(), bound_insert_location(left)); +#ifdef GCC_MISSING_VECTOR_RANGE_INSERT + itr = active_bounds.insert(itr, &right); + return active_bounds.insert(itr, &left); +#else + return active_bounds.insert(itr, { &left, &right }); +#endif +} + +template +inline bool is_maxima(bound const& bnd, T y) { + return bnd.next_edge == bnd.edges.end() && bnd.current_edge->top.y == y; +} + +template +inline bool is_maxima(active_bound_list_itr const& bnd, T y) { + return is_maxima(*(*bnd), y); +} + +template +inline bool is_intermediate(bound const& bnd, T y) { + return bnd.next_edge != bnd.edges.end() && bnd.current_edge->top.y == y; +} + +template +inline bool is_intermediate(active_bound_list_itr const& bnd, T y) { + return is_intermediate(*(*bnd), y); +} + +template +inline bool current_edge_is_horizontal(active_bound_list_itr const& bnd) { + return is_horizontal(*((*bnd)->current_edge)); +} + +template +inline bool next_edge_is_horizontal(active_bound_list_itr const& bnd) { + return is_horizontal(*((*bnd)->next_edge)); +} + +template +void next_edge_in_bound(bound& bnd, scanbeam_list& scanbeam) { + auto& current_edge = bnd.current_edge; + ++current_edge; + if (current_edge != bnd.edges.end()) { + ++(bnd.next_edge); + bnd.current_x = static_cast(current_edge->bot.x); + if (!is_horizontal(*current_edge)) { + insert_sorted_scanbeam(scanbeam, current_edge->top.y); + } + } +} + +template +active_bound_list_itr get_maxima_pair(active_bound_list_itr bnd, active_bound_list& active_bounds) { + bound_ptr maximum = (*bnd)->maximum_bound; + return std::find(active_bounds.begin(), active_bounds.end(), maximum); +} + +template +void set_winding_count(active_bound_list_itr bnd_itr, + active_bound_list& active_bounds, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + auto rev_bnd_itr = active_bound_list_rev_itr(bnd_itr); + if (rev_bnd_itr == active_bounds.rend()) { + (*bnd_itr)->winding_count = (*bnd_itr)->winding_delta; + (*bnd_itr)->winding_count2 = 0; + return; + } + + // find the edge of the same polytype that immediately preceeds 'edge' in + // AEL + while (rev_bnd_itr != active_bounds.rend() && (*rev_bnd_itr)->poly_type != (*bnd_itr)->poly_type) { + ++rev_bnd_itr; + } + if (rev_bnd_itr == active_bounds.rend()) { + (*bnd_itr)->winding_count = (*bnd_itr)->winding_delta; + (*bnd_itr)->winding_count2 = 0; + } else if (is_even_odd_fill_type(*(*bnd_itr), subject_fill_type, clip_fill_type)) { + // EvenOdd filling ... + (*bnd_itr)->winding_count = (*bnd_itr)->winding_delta; + (*bnd_itr)->winding_count2 = (*rev_bnd_itr)->winding_count2; + } else { + // nonZero, Positive or Negative filling ... + if ((*rev_bnd_itr)->winding_count * (*rev_bnd_itr)->winding_delta < 0) { + // prev edge is 'decreasing' WindCount (WC) toward zero + // so we're outside the previous polygon ... + if (std::abs(static_cast((*rev_bnd_itr)->winding_count)) > 1) { + // outside prev poly but still inside another. + // when reversing direction of prev poly use the same WC + if ((*rev_bnd_itr)->winding_delta * (*bnd_itr)->winding_delta < 0) { + (*bnd_itr)->winding_count = (*rev_bnd_itr)->winding_count; + } else { + // otherwise continue to 'decrease' WC ... + (*bnd_itr)->winding_count = (*rev_bnd_itr)->winding_count + (*bnd_itr)->winding_delta; + } + } else { + // now outside all polys of same polytype so set own WC ... + (*bnd_itr)->winding_count = (*bnd_itr)->winding_delta; + } + } else { + // prev edge is 'increasing' WindCount (WC) away from zero + // so we're inside the previous polygon ... + if ((*rev_bnd_itr)->winding_delta * (*bnd_itr)->winding_delta < 0) { + // if wind direction is reversing prev then use same WC + (*bnd_itr)->winding_count = (*rev_bnd_itr)->winding_count; + } else { + // otherwise add to WC ... + (*bnd_itr)->winding_count = (*rev_bnd_itr)->winding_count + (*bnd_itr)->winding_delta; + } + } + (*bnd_itr)->winding_count2 = (*rev_bnd_itr)->winding_count2; + } + + // update winding_count2 ... + auto bnd_itr_forward = rev_bnd_itr.base(); + if (is_even_odd_alt_fill_type(*(*bnd_itr), subject_fill_type, clip_fill_type)) { + // EvenOdd filling ... + while (bnd_itr_forward != bnd_itr) { + (*bnd_itr)->winding_count2 = ((*bnd_itr)->winding_count2 == 0 ? 1 : 0); + ++bnd_itr_forward; + } + } else { + // nonZero, Positive or Negative filling ... + while (bnd_itr_forward != bnd_itr) { + (*bnd_itr)->winding_count2 += (*bnd_itr_forward)->winding_delta; + ++bnd_itr_forward; + } + } +} + +template +bool is_contributing(bound const& bnd, clip_type cliptype, fill_type subject_fill_type, fill_type clip_fill_type) { + fill_type pft = subject_fill_type; + fill_type pft2 = clip_fill_type; + if (bnd.poly_type != polygon_type_subject) { + pft = clip_fill_type; + pft2 = subject_fill_type; + } + + switch (pft) { + case fill_type_even_odd: + break; + case fill_type_non_zero: + if (std::abs(static_cast(bnd.winding_count)) != 1) { + return false; + } + break; + case fill_type_positive: + if (bnd.winding_count != 1) { + return false; + } + break; + case fill_type_negative: + default: + if (bnd.winding_count != -1) { + return false; + } + } + + switch (cliptype) { + case clip_type_intersection: + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 != 0); + case fill_type_positive: + return (bnd.winding_count2 > 0); + case fill_type_negative: + default: + return (bnd.winding_count2 < 0); + } + break; + case clip_type_union: + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 == 0); + case fill_type_positive: + return (bnd.winding_count2 <= 0); + case fill_type_negative: + default: + return (bnd.winding_count2 >= 0); + } + break; + case clip_type_difference: + if (bnd.poly_type == polygon_type_subject) { + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 == 0); + case fill_type_positive: + return (bnd.winding_count2 <= 0); + case fill_type_negative: + default: + return (bnd.winding_count2 >= 0); + } + } else { + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 != 0); + case fill_type_positive: + return (bnd.winding_count2 > 0); + case fill_type_negative: + default: + return (bnd.winding_count2 < 0); + } + } + break; + case clip_type_x_or: + return true; + break; + default: + return true; + } +} + +template +void insert_lm_left_and_right_bound(bound& left_bound, + bound& right_bound, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + // Both left and right bound + auto lb_abl_itr = insert_bound_into_ABL(left_bound, right_bound, active_bounds); + auto rb_abl_itr = std::next(lb_abl_itr); + set_winding_count(lb_abl_itr, active_bounds, subject_fill_type, clip_fill_type); + (*rb_abl_itr)->winding_count = (*lb_abl_itr)->winding_count; + (*rb_abl_itr)->winding_count2 = (*lb_abl_itr)->winding_count2; + if (is_contributing(left_bound, cliptype, subject_fill_type, clip_fill_type)) { + add_local_minimum_point(*(*lb_abl_itr), *(*rb_abl_itr), active_bounds, (*lb_abl_itr)->current_edge->bot, rings); + } + + // Add top of edges to scanbeam + insert_sorted_scanbeam(scanbeam, (*lb_abl_itr)->current_edge->top.y); + + if (!current_edge_is_horizontal(rb_abl_itr)) { + insert_sorted_scanbeam(scanbeam, (*rb_abl_itr)->current_edge->top.y); + } +} + +template +void insert_local_minima_into_ABL(T const bot_y, + local_minimum_ptr_list const& minima_sorted, + local_minimum_ptr_list_itr& current_lm, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + while (current_lm != minima_sorted.end() && bot_y == (*current_lm)->y) { + initialize_lm(current_lm); + auto& left_bound = (*current_lm)->left_bound; + auto& right_bound = (*current_lm)->right_bound; + insert_lm_left_and_right_bound(left_bound, right_bound, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + ++current_lm; + } +} + +template +void insert_horizontal_local_minima_into_ABL(T const top_y, + local_minimum_ptr_list const& minima_sorted, + local_minimum_ptr_list_itr& current_lm, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + while (current_lm != minima_sorted.end() && top_y == (*current_lm)->y && (*current_lm)->minimum_has_horizontal) { + initialize_lm(current_lm); + auto& left_bound = (*current_lm)->left_bound; + auto& right_bound = (*current_lm)->right_bound; + insert_lm_left_and_right_bound(left_bound, right_bound, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + ++current_lm; + } +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/bound.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/bound.hpp new file mode 100644 index 000000000..fa748fdaf --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/bound.hpp @@ -0,0 +1,99 @@ +#pragma once + +#include + +#include + +#include +#include +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct bound { + + edge_list edges; + edge_list_itr current_edge; + edge_list_itr next_edge; + mapbox::geometry::point last_point; + ring_ptr ring; + bound_ptr maximum_bound; // the bound who's maximum connects with this bound + double current_x; + std::size_t pos; + std::int32_t winding_count; + std::int32_t winding_count2; // winding count of the opposite polytype + std::int8_t winding_delta; // 1 or -1 depending on winding direction - 0 for linestrings + polygon_type poly_type; + edge_side side; // side only refers to current side of solution poly + + bound() noexcept + : edges(), + current_edge(edges.end()), + next_edge(edges.end()), + last_point({ 0, 0 }), + ring(nullptr), + maximum_bound(nullptr), + current_x(0.0), + pos(0), + winding_count(0), + winding_count2(0), + winding_delta(0), + poly_type(polygon_type_subject), + side(edge_left) { + } + + bound(bound&& b) noexcept + : edges(std::move(b.edges)), + current_edge(std::move(b.current_edge)), + next_edge(std::move(b.next_edge)), + last_point(std::move(b.last_point)), + ring(std::move(b.ring)), + maximum_bound(std::move(b.maximum_bound)), + current_x(std::move(b.current_x)), + pos(std::move(b.pos)), + winding_count(std::move(b.winding_count)), + winding_count2(std::move(b.winding_count2)), + winding_delta(std::move(b.winding_delta)), + poly_type(std::move(b.poly_type)), + side(std::move(b.side)) { + } + + bound(bound const& b) = delete; + bound& operator=(bound const&) = delete; +}; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, const bound& bnd) { + out << " Bound: " << &bnd << std::endl; + out << " current_x: " << bnd.current_x << std::endl; + out << " last_point: " << bnd.last_point.x << ", " << bnd.last_point.y << std::endl; + out << *(bnd.current_edge); + out << " winding count: " << bnd.winding_count << std::endl; + out << " winding_count2: " << bnd.winding_count2 << std::endl; + out << " winding_delta: " << static_cast(bnd.winding_delta) << std::endl; + out << " maximum_bound: " << bnd.maximum_bound << std::endl; + if (bnd.side == edge_left) { + out << " side: left" << std::endl; + } else { + out << " side: right" << std::endl; + } + out << " ring: " << bnd.ring << std::endl; + if (bnd.ring) { + out << " ring index: " << bnd.ring->ring_index << std::endl; + } + return out; +} + +#endif +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/bubble_sort.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/bubble_sort.hpp new file mode 100644 index 000000000..75f8d212b --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/bubble_sort.hpp @@ -0,0 +1,28 @@ +#pragma once + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +void bubble_sort(It begin, It end, Compare c, MethodOnSwap m) { + if (begin == end) { + return; + } + bool modified = false; + auto last = end - 1; + do { + modified = false; + for (auto itr = begin; itr != last; ++itr) { + auto next = std::next(itr); + if (!c(*itr, *next)) { + m(*itr, *next); + std::iter_swap(itr, next); + modified = true; + } + } + } while (modified); +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/build_edges.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/build_edges.hpp new file mode 100644 index 000000000..a1131d8be --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/build_edges.hpp @@ -0,0 +1,183 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +bool point_2_is_between_point_1_and_point_3(mapbox::geometry::point const& pt1, + mapbox::geometry::point const& pt2, + mapbox::geometry::point const& pt3) { + if ((pt1 == pt3) || (pt1 == pt2) || (pt3 == pt2)) { + return false; + } else if (pt1.x != pt3.x) { + return (pt2.x > pt1.x) == (pt2.x < pt3.x); + } else { + return (pt2.y > pt1.y) == (pt2.y < pt3.y); + } +} + +template +bool build_edge_list(mapbox::geometry::linear_ring const& path_geometry, edge_list& edges) { + + if (path_geometry.size() < 3) { + return false; + } + + // As this is a loop, we need to first go backwards from end to try and find + // the proper starting point for the iterators before the beginning + + auto itr_rev = path_geometry.rbegin(); + auto itr = path_geometry.begin(); + mapbox::geometry::point pt1 = *itr_rev; + mapbox::geometry::point pt2 = *itr; + + // Find next non repeated point going backwards from + // end for pt1 + while (pt1 == pt2) { + ++itr_rev; + if (itr_rev == path_geometry.rend()) { + return false; + } + pt1 = *itr_rev; + } + ++itr; + mapbox::geometry::point pt3 = *itr; + auto itr_last = itr_rev.base(); + mapbox::geometry::point front_pt; + mapbox::geometry::point back_pt; + while (true) { + if (pt3 == pt2) { + // Duplicate point advance itr, but do not + // advance other points + if (itr == itr_last) { + break; + } + ++itr; + if (itr == itr_last) { + if (edges.empty()) { + break; + } + pt3 = front_pt; + } else { + pt3 = *itr; + } + continue; + } + + // Now check if slopes are equal between two segments - either + // a spike or a collinear point - if so drop point number 2. + if (slopes_equal(pt1, pt2, pt3)) { + // We need to reconsider previously added points + // because the point it was using was found to be collinear + // or a spike + pt2 = pt1; + if (!edges.empty()) { + edges.pop_back(); // remove previous edge (pt1) + } + if (!edges.empty()) { + auto const& back_top = edges.back().top; + if (static_cast(back_pt.x) == back_top.x && static_cast(back_pt.y) == back_top.y) { + auto const& back_bot = edges.back().bot; + pt1 = mapbox::geometry::point(static_cast(back_bot.x), static_cast(back_bot.y)); + } else { + pt1 = mapbox::geometry::point(static_cast(back_top.x), static_cast(back_top.y)); + } + back_pt = pt1; + } else { + // If this occurs we must look to the back of the + // ring for new points. + while (*itr_rev == pt2) { + ++itr_rev; + if ((itr + 1) == itr_rev.base()) { + return false; + } + } + pt1 = *itr_rev; + itr_last = itr_rev.base(); + } + continue; + } + + if (edges.empty()) { + front_pt = pt2; + } + edges.emplace_back(pt2, pt3); + back_pt = pt2; + if (itr == itr_last) { + break; + } + pt1 = pt2; + pt2 = pt3; + ++itr; + if (itr == itr_last) { + if (edges.empty()) { + break; + } + pt3 = front_pt; + } else { + pt3 = *itr; + } + } + + bool modified = false; + do { + modified = false; + if (edges.size() < 3) { + return false; + } + auto& f = edges.front(); + auto& b = edges.back(); + if (slopes_equal(f, b)) { + if (f.bot == b.top) { + if (f.top == b.bot) { + edges.pop_back(); + edges.erase(edges.begin()); + } else { + f.bot = b.bot; + edges.pop_back(); + } + modified = true; + } else if (f.top == b.bot) { + f.top = b.top; + edges.pop_back(); + modified = true; + } else if (f.top == b.top && f.bot == b.bot) { + edges.pop_back(); + edges.erase(edges.begin()); + modified = true; + } else if (f.top == b.top) { + if (point_2_is_between_point_1_and_point_3(f.top, f.bot, b.bot)) { + b.top = f.bot; + edges.erase(edges.begin()); + } else { + f.top = b.bot; + edges.pop_back(); + } + modified = true; + } else if (f.bot == b.bot) { + if (point_2_is_between_point_1_and_point_3(f.bot, f.top, b.top)) { + b.bot = f.top; + edges.erase(edges.begin()); + } else { + f.bot = b.top; + edges.pop_back(); + } + modified = true; + } + } + } while (modified); + + return true; +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/build_local_minima_list.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/build_local_minima_list.hpp new file mode 100644 index 000000000..3ed42a5b1 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/build_local_minima_list.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +bool add_linear_ring(mapbox::geometry::linear_ring const& path_geometry, + local_minimum_list& minima_list, + polygon_type p_type) { + edge_list new_edges; + new_edges.reserve(path_geometry.size()); + if (!build_edge_list(path_geometry, new_edges) || new_edges.empty()) { + return false; + } + add_ring_to_local_minima_list(new_edges, minima_list, p_type); + return true; +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/build_result.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/build_result.hpp new file mode 100644 index 000000000..261167850 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/build_result.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include +#include + +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +void push_ring_to_polygon(mapbox::geometry::polygon& poly, ring_ptr r, bool reverse_output) { + mapbox::geometry::linear_ring lr; + lr.reserve(r->size() + 1); + auto firstPt = r->points; + auto ptIt = r->points; + if (reverse_output) { + do { + lr.emplace_back(static_cast(ptIt->x), static_cast(ptIt->y)); + ptIt = ptIt->next; + } while (ptIt != firstPt); + } else { + do { + lr.emplace_back(static_cast(ptIt->x), static_cast(ptIt->y)); + ptIt = ptIt->prev; + } while (ptIt != firstPt); + } + lr.emplace_back(firstPt->x, firstPt->y); // close the ring + poly.push_back(lr); +} + +template +void build_result_polygons(mapbox::geometry::multi_polygon& solution, + ring_vector const& rings, + bool reverse_output) { + for (auto r : rings) { + if (r == nullptr) { + continue; + } + assert(r->points); + solution.emplace_back(); + push_ring_to_polygon(solution.back(), r, reverse_output); + for (auto c : r->children) { + if (c == nullptr) { + continue; + } + assert(c->points); + push_ring_to_polygon(solution.back(), c, reverse_output); + } + for (auto c : r->children) { + if (c == nullptr) { + continue; + } + if (!c->children.empty()) { + build_result_polygons(solution, c->children, reverse_output); + } + } + } +} + +template +void build_result(mapbox::geometry::multi_polygon& solution, ring_manager const& rings, bool reverse_output) { + build_result_polygons(solution, rings.children, reverse_output); +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/config.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/config.hpp new file mode 100644 index 000000000..1d2211a26 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/config.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include +#include +#include +#include + +// GCC 4.8 missing range std::vector::insert (c++11) +#ifdef __GNUC__ +#if __GNUC__ == 4 && __GNUC_MINOR__ == 8 +#define GCC_MISSING_VECTOR_RANGE_INSERT +#endif +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +enum clip_type : std::uint8_t { clip_type_intersection = 0, clip_type_union, clip_type_difference, clip_type_x_or }; + +enum polygon_type : std::uint8_t { polygon_type_subject = 0, polygon_type_clip }; + +enum fill_type : std::uint8_t { fill_type_even_odd = 0, fill_type_non_zero, fill_type_positive, fill_type_negative }; + +static double const def_arc_tolerance = 0.25; + +static int const EDGE_UNASSIGNED = -1; // edge not currently 'owning' a solution +static int const EDGE_SKIP = -2; // edge that would otherwise close a path +static std::int64_t const LOW_RANGE = 0x3FFFFFFF; +static std::int64_t const HIGH_RANGE = 0x3FFFFFFFFFFFFFFFLL; + +enum horizontal_direction : std::uint8_t { right_to_left = 0, left_to_right = 1 }; + +enum edge_side : std::uint8_t { edge_left = 0, edge_right }; + +enum join_type : std::uint8_t { join_type_square = 0, join_type_round, join_type_miter }; + +enum end_type { + end_type_closed_polygon = 0, + end_type_closed_line, + end_type_open_butt, + end_type_open_square, + end_type_open_round +}; + +template +using maxima_list = std::list; +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/edge.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/edge.hpp new file mode 100644 index 000000000..6178cb660 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/edge.hpp @@ -0,0 +1,120 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct bound; + +template +using bound_ptr = bound*; + +template +struct edge { + mapbox::geometry::point bot; + mapbox::geometry::point top; + double dx; + + edge(edge&& e) noexcept : bot(std::move(e.bot)), top(std::move(e.top)), dx(std::move(e.dx)) { + } + + edge& operator=(edge&& e) noexcept { + bot = std::move(e.bot); + top = std::move(e.top); + dx = std::move(e.dx); + return *this; + } + + template + edge(mapbox::geometry::point const& current, mapbox::geometry::point const& next_pt) noexcept + : bot(static_cast(current.x), static_cast(current.y)), + top(static_cast(current.x), static_cast(current.y)), + dx(0.0) { + if (current.y >= next_pt.y) { + top = mapbox::geometry::point(static_cast(next_pt.x), static_cast(next_pt.y)); + } else { + bot = mapbox::geometry::point(static_cast(next_pt.x), static_cast(next_pt.y)); + } + double dy = static_cast(top.y - bot.y); + if (value_is_zero(dy)) { + dx = std::numeric_limits::infinity(); + } else { + dx = static_cast(top.x - bot.x) / dy; + } + } +}; + +template +using edge_ptr = edge*; + +template +using edge_list = std::vector>; + +template +using edge_list_itr = typename edge_list::iterator; + +template +bool slopes_equal(edge const& e1, edge const& e2) { + return static_cast(e1.top.y - e1.bot.y) * static_cast(e2.top.x - e2.bot.x) == + static_cast(e1.top.x - e1.bot.x) * static_cast(e2.top.y - e2.bot.y); +} + +template +inline bool is_horizontal(edge const& e) { + return std::isinf(e.dx); +} + +template +inline double get_current_x(edge const& edge, const T current_y) { + if (current_y == edge.top.y) { + return static_cast(edge.top.x); + } else { + return static_cast(edge.bot.x) + edge.dx * static_cast(current_y - edge.bot.y); + } +} + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, const edge& e) { + out << " Edge: " << std::endl; + out << " bot x: " << e.bot.x << " y: " << e.bot.y << std::endl; + out << " top x: " << e.top.x << " y: " << e.top.y << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + edge_list const& edges) { + out << "["; + bool first = true; + for (auto const& e : edges) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << e.bot.x << "," << e.bot.y << "],["; + out << e.top.x << "," << e.top.y << "]]"; + } + out << "]"; + return out; +} + +#endif +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/interrupt.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/interrupt.hpp new file mode 100644 index 000000000..083cf360d --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/interrupt.hpp @@ -0,0 +1,50 @@ +#pragma once + +/** + * To enable this by the program, define USE_WAGYU_INTERRUPT before including wagyu.hpp + * To request an interruption, call `interrupt_request()`. As soon as Wagyu detects the request + * it will raise an exception (`std::runtime_error`). + */ + +#ifdef USE_WAGYU_INTERRUPT + +namespace { +bool WAGYU_INTERRUPT_REQUESTED = false; +} + +namespace mapbox { +namespace geometry { +namespace wagyu { + +static void interrupt_reset(void) { + WAGYU_INTERRUPT_REQUESTED = false; +} + +static void interrupt_request(void) { + WAGYU_INTERRUPT_REQUESTED = true; +} + +static void interrupt_check(void) { + if (WAGYU_INTERRUPT_REQUESTED) { + interrupt_reset(); + throw std::runtime_error("Wagyu interrupted"); + } +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox + +#else /* ! USE_WAGYU_INTERRUPT */ + +namespace mapbox { +namespace geometry { +namespace wagyu { + +static void interrupt_check(void) { +} + +} // namespace wagyu +} // namespace geometry +} // namespace mapbox + +#endif /* USE_WAGYU_INTERRUPT */ diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/intersect.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/intersect.hpp new file mode 100644 index 000000000..aae4c8f98 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/intersect.hpp @@ -0,0 +1,70 @@ +#pragma once + +#include + +#include + +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct intersect_node { + + bound_ptr bound1; + bound_ptr bound2; + mapbox::geometry::point pt; + + intersect_node(intersect_node&& n) noexcept + : bound1(std::move(n.bound1)), bound2(std::move(n.bound2)), pt(std::move(n.pt)) { + } + + intersect_node& operator=(intersect_node&& n) noexcept { + bound1 = std::move(n.bound1); + bound2 = std::move(n.bound2); + pt = std::move(n.pt); + return *this; + } + + intersect_node(bound_ptr const& bound1_, bound_ptr const& bound2_, mapbox::geometry::point const& pt_) + : bound1(bound1_), bound2(bound2_), pt(pt_) { + } +}; + +template +using intersect_list = std::vector>; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const intersect_node& e) { + out << " point x: " << e.pt.x << " y: " << e.pt.y << std::endl; + out << " bound 1: " << std::endl; + out << *e.bound1 << std::endl; + out << " bound 2: " << std::endl; + out << *e.bound2 << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const intersect_list& ints) { + std::size_t c = 0; + for (auto const& i : ints) { + out << "Intersection: " << c++ << std::endl; + out << i; + } + return out; +} + +#endif +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/intersect_util.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/intersect_util.hpp new file mode 100644 index 000000000..34dd35497 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/intersect_util.hpp @@ -0,0 +1,365 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct intersect_list_sorter { + inline bool operator()(intersect_node const& node1, intersect_node const& node2) { + if (!values_are_equal(node2.pt.y, node1.pt.y)) { + return node2.pt.y < node1.pt.y; + } else { + return (node2.bound1->winding_count2 + node2.bound2->winding_count2) > + (node1.bound1->winding_count2 + node1.bound2->winding_count2); + } + } +}; + +template +inline mapbox::geometry::point round_point(mapbox::geometry::point const& pt) { + return mapbox::geometry::point(round_towards_max(pt.x), round_towards_max(pt.y)); +} + +template +inline void swap_rings(bound& b1, bound& b2) { + ring_ptr ring = b1.ring; + b1.ring = b2.ring; + b2.ring = ring; +} + +template +inline void swap_sides(bound& b1, bound& b2) { + edge_side side = b1.side; + b1.side = b2.side; + b2.side = side; +} + +template +bool get_edge_intersection(edge const& e1, edge const& e2, mapbox::geometry::point& pt) { + T2 p0_x = static_cast(e1.bot.x); + T2 p0_y = static_cast(e1.bot.y); + T2 p1_x = static_cast(e1.top.x); + T2 p1_y = static_cast(e1.top.y); + T2 p2_x = static_cast(e2.bot.x); + T2 p2_y = static_cast(e2.bot.y); + T2 p3_x = static_cast(e2.top.x); + T2 p3_y = static_cast(e2.top.y); + T2 s1_x, s1_y, s2_x, s2_y; + s1_x = p1_x - p0_x; + s1_y = p1_y - p0_y; + s2_x = p3_x - p2_x; + s2_y = p3_y - p2_y; + + T2 s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / (-s2_x * s1_y + s1_x * s2_y); + T2 t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / (-s2_x * s1_y + s1_x * s2_y); + + if (s >= 0.0 && s <= 1.0 && t >= 0.0 && t <= 1.0) { + pt.x = p0_x + (t * s1_x); + pt.y = p0_y + (t * s1_y); + return true; + } + // LCOV_EXCL_START + return false; + // LCOV_EXCL_END +} + +template +struct intersection_compare { + bool operator()(bound_ptr const& b1, bound_ptr const& b2) { + return !(b1->current_x > b2->current_x && !slopes_equal(*(b1->current_edge), *(b2->current_edge))); + } +}; + +template +struct on_intersection_swap { + + intersect_list& intersects; + + on_intersection_swap(intersect_list& i) : intersects(i) { + } + + void operator()(bound_ptr const& b1, bound_ptr const& b2) { + mapbox::geometry::point pt; + if (!get_edge_intersection(*(b1->current_edge), *(b2->current_edge), pt)) { + // LCOV_EXCL_START + throw std::runtime_error("Trying to find intersection of lines that do not intersect"); + // LCOV_EXCL_END + } + intersects.emplace_back(b1, b2, pt); + } +}; + +template +void build_intersect_list(active_bound_list& active_bounds, intersect_list& intersects) { + bubble_sort(active_bounds.begin(), active_bounds.end(), intersection_compare(), + on_intersection_swap(intersects)); +} + +template +void intersect_bounds(bound& b1, + bound& b2, + mapbox::geometry::point const& pt, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& rings, + active_bound_list& active_bounds) { + bool b1Contributing = (b1.ring != nullptr); + bool b2Contributing = (b2.ring != nullptr); + + // update winding counts... + // assumes that b1 will be to the Right of b2 ABOVE the intersection + if (b1.poly_type == b2.poly_type) { + if (is_even_odd_fill_type(b1, subject_fill_type, clip_fill_type)) { + std::swap(b1.winding_count, b2.winding_count); + } else { + if (b1.winding_count + b2.winding_delta == 0) { + b1.winding_count = -b1.winding_count; + } else { + b1.winding_count += b2.winding_delta; + } + if (b2.winding_count - b1.winding_delta == 0) { + b2.winding_count = -b2.winding_count; + } else { + b2.winding_count -= b1.winding_delta; + } + } + } else { + if (!is_even_odd_fill_type(b2, subject_fill_type, clip_fill_type)) { + b1.winding_count2 += b2.winding_delta; + } else { + b1.winding_count2 = (b1.winding_count2 == 0) ? 1 : 0; + } + if (!is_even_odd_fill_type(b1, subject_fill_type, clip_fill_type)) { + b2.winding_count2 -= b1.winding_delta; + } else { + b2.winding_count2 = (b2.winding_count2 == 0) ? 1 : 0; + } + } + + fill_type b1FillType, b2FillType, b1FillType2, b2FillType2; + if (b1.poly_type == polygon_type_subject) { + b1FillType = subject_fill_type; + b1FillType2 = clip_fill_type; + } else { + b1FillType = clip_fill_type; + b1FillType2 = subject_fill_type; + } + if (b2.poly_type == polygon_type_subject) { + b2FillType = subject_fill_type; + b2FillType2 = clip_fill_type; + } else { + b2FillType = clip_fill_type; + b2FillType2 = subject_fill_type; + } + + std::int32_t b1Wc, b2Wc; + switch (b1FillType) { + case fill_type_positive: + b1Wc = b1.winding_count; + break; + case fill_type_negative: + b1Wc = -b1.winding_count; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b1Wc = std::abs(static_cast(b1.winding_count)); + } + switch (b2FillType) { + case fill_type_positive: + b2Wc = b2.winding_count; + break; + case fill_type_negative: + b2Wc = -b2.winding_count; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b2Wc = std::abs(static_cast(b2.winding_count)); + } + if (b1Contributing && b2Contributing) { + if ((b1Wc != 0 && b1Wc != 1) || (b2Wc != 0 && b2Wc != 1) || + (b1.poly_type != b2.poly_type && cliptype != clip_type_x_or)) { + add_local_maximum_point(b1, b2, pt, rings, active_bounds); + } else { + add_point(b1, active_bounds, pt, rings); + add_point(b2, active_bounds, pt, rings); + swap_sides(b1, b2); + swap_rings(b1, b2); + } + } else if (b1Contributing) { + if (b2Wc == 0 || b2Wc == 1) { + add_point(b1, active_bounds, pt, rings); + b2.last_point = pt; + swap_sides(b1, b2); + swap_rings(b1, b2); + } + } else if (b2Contributing) { + if (b1Wc == 0 || b1Wc == 1) { + b1.last_point = pt; + add_point(b2, active_bounds, pt, rings); + swap_sides(b1, b2); + swap_rings(b1, b2); + } + } else if ((b1Wc == 0 || b1Wc == 1) && (b2Wc == 0 || b2Wc == 1)) { + // neither bound is currently contributing ... + + std::int32_t b1Wc2, b2Wc2; + switch (b1FillType2) { + case fill_type_positive: + b1Wc2 = b1.winding_count2; + break; + case fill_type_negative: + b1Wc2 = -b1.winding_count2; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b1Wc2 = std::abs(static_cast(b1.winding_count2)); + } + switch (b2FillType2) { + case fill_type_positive: + b2Wc2 = b2.winding_count2; + break; + case fill_type_negative: + b2Wc2 = -b2.winding_count2; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b2Wc2 = std::abs(static_cast(b2.winding_count2)); + } + + if (b1.poly_type != b2.poly_type) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } else if (b1Wc == 1 && b2Wc == 1) { + switch (cliptype) { + case clip_type_intersection: + if (b1Wc2 > 0 && b2Wc2 > 0) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + break; + default: + case clip_type_union: + if (b1Wc2 <= 0 && b2Wc2 <= 0) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + break; + case clip_type_difference: + if (((b1.poly_type == polygon_type_clip) && (b1Wc2 > 0) && (b2Wc2 > 0)) || + ((b1.poly_type == polygon_type_subject) && (b1Wc2 <= 0) && (b2Wc2 <= 0))) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + break; + case clip_type_x_or: + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + } else { + swap_sides(b1, b2); + } + } +} + +template +bool bounds_adjacent(intersect_node const& inode, bound_ptr next) { + return (next == inode.bound2) || (next == inode.bound1); +} + +template +struct find_first_bound { + bound_ptr b1; + bound_ptr b2; + + find_first_bound(intersect_node const& inode) : b1(inode.bound1), b2(inode.bound2) { + } + + bool operator()(bound_ptr const& b) { + return b == b1 || b == b2; + } +}; + +template +void process_intersect_list(intersect_list& intersects, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& rings, + active_bound_list& active_bounds) { + for (auto node_itr = intersects.begin(); node_itr != intersects.end(); ++node_itr) { + auto b1 = std::find_if(active_bounds.begin(), active_bounds.end(), find_first_bound(*node_itr)); + auto b2 = std::next(b1); + if (!bounds_adjacent(*node_itr, *b2)) { + auto next_itr = std::next(node_itr); + while (next_itr != intersects.end()) { + auto n1 = std::find_if(active_bounds.begin(), active_bounds.end(), find_first_bound(*next_itr)); + auto n2 = std::next(n1); + if (bounds_adjacent(*next_itr, *n2)) { + b1 = n1; + b2 = n2; + break; + } + ++next_itr; + } + if (next_itr == intersects.end()) { + throw std::runtime_error("Could not properly correct intersection order."); + } + std::iter_swap(node_itr, next_itr); + } + mapbox::geometry::point pt = round_point(node_itr->pt); + intersect_bounds(*(node_itr->bound1), *(node_itr->bound2), pt, cliptype, subject_fill_type, clip_fill_type, + rings, active_bounds); + std::iter_swap(b1, b2); + } +} + +template +void update_current_x(active_bound_list& active_bounds, T top_y) { + std::size_t pos = 0; + for (auto& bnd : active_bounds) { + bnd->pos = pos++; + bnd->current_x = get_current_x(*bnd->current_edge, top_y); + } +} + +template +void process_intersections(T top_y, + active_bound_list& active_bounds, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& rings) { + if (active_bounds.empty()) { + return; + } + update_current_x(active_bounds, top_y); + intersect_list intersects; + build_intersect_list(active_bounds, intersects); + + if (intersects.empty()) { + return; + } + + // Restore order of active bounds list + std::stable_sort(active_bounds.begin(), active_bounds.end(), + [](bound_ptr const& b1, bound_ptr const& b2) { return b1->pos < b2->pos; }); + + // Sort the intersection list + std::stable_sort(intersects.begin(), intersects.end(), intersect_list_sorter()); + + process_intersect_list(intersects, cliptype, subject_fill_type, clip_fill_type, rings, active_bounds); +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/local_minimum.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/local_minimum.hpp new file mode 100644 index 000000000..1f579106d --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/local_minimum.hpp @@ -0,0 +1,117 @@ +#pragma once + +#ifdef DEBUG +#include +#include +#endif + +#include + +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct local_minimum { + bound left_bound; + bound right_bound; + T y; + bool minimum_has_horizontal; + + local_minimum(bound&& left_bound_, bound&& right_bound_, T y_, bool has_horz_) + : left_bound(std::move(left_bound_)), + right_bound(std::move(right_bound_)), + y(y_), + minimum_has_horizontal(has_horz_) { + } +}; + +template +using local_minimum_list = std::deque>; + +template +using local_minimum_itr = typename local_minimum_list::iterator; + +template +using local_minimum_ptr = local_minimum*; + +template +using local_minimum_ptr_list = std::vector>; + +template +using local_minimum_ptr_list_itr = typename local_minimum_ptr_list::iterator; + +template +struct local_minimum_sorter { + inline bool operator()(local_minimum_ptr const& locMin1, local_minimum_ptr const& locMin2) { + if (locMin2->y == locMin1->y) { + return locMin2->minimum_has_horizontal != locMin1->minimum_has_horizontal && + locMin1->minimum_has_horizontal; + } + return locMin2->y < locMin1->y; + } +}; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const local_minimum& lm) { + out << " Local Minimum:" << std::endl; + out << " y: " << lm.y << std::endl; + if (lm.minimum_has_horizontal) { + out << " minimum_has_horizontal: true" << std::endl; + } else { + out << " minimum_has_horizontal: false" << std::endl; + } + out << " left_bound: " << std::endl; + out << lm.left_bound << std::endl; + out << " right_bound: " << std::endl; + out << lm.right_bound << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const local_minimum_ptr_list& lms) { + for (auto const& lm : lms) { + out << *lm; + } + return out; +} + +template +std::string output_all_edges(local_minimum_ptr_list const& lms) { + std::ostringstream out; + out << "["; + bool first = true; + for (auto const& lm : lms) { + for (auto const& e : lm->left_bound.edges) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << e.bot.x << "," << e.bot.y << "],["; + out << e.top.x << "," << e.top.y << "]]"; + } + for (auto const& e : lm->right_bound.edges) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << e.bot.x << "," << e.bot.y << "],["; + out << e.top.x << "," << e.top.y << "]]"; + } + } + out << "]"; + return out.str(); +} + +#endif +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/local_minimum_util.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/local_minimum_util.hpp new file mode 100644 index 000000000..50fd0f563 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/local_minimum_util.hpp @@ -0,0 +1,314 @@ +#pragma once + +#include +#include +#include + +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +inline void reverse_horizontal(edge& e) { + // swap horizontal edges' top and bottom x's so they follow the natural + // progression of the bounds - ie so their xbots will align with the + // adjoining lower edge. [Helpful in the process_horizontal() method.] + std::swap(e.top.x, e.bot.x); +} + +// Make a list start on a local maximum by +// shifting all the points not on a local maximum to the +template +void start_list_on_local_maximum(edge_list& edges) { + if (edges.size() <= 2) { + return; + } + // Find the first local maximum going forward in the list + auto prev_edge = edges.end(); + --prev_edge; + bool prev_edge_is_horizontal = is_horizontal(*prev_edge); + auto edge = edges.begin(); + bool edge_is_horizontal; + bool y_decreasing_before_last_horizontal = false; // assume false at start + + while (edge != edges.end()) { + edge_is_horizontal = is_horizontal(*edge); + if ((!prev_edge_is_horizontal && !edge_is_horizontal && edge->top == prev_edge->top)) { + break; + } + if (!edge_is_horizontal && prev_edge_is_horizontal) { + if (y_decreasing_before_last_horizontal && (edge->top == prev_edge->bot || edge->top == prev_edge->top)) { + break; + } + } else if (!y_decreasing_before_last_horizontal && !prev_edge_is_horizontal && edge_is_horizontal && + (prev_edge->top == edge->top || prev_edge->top == edge->bot)) { + y_decreasing_before_last_horizontal = true; + } + prev_edge_is_horizontal = edge_is_horizontal; + prev_edge = edge; + ++edge; + } + std::rotate(edges.begin(), edge, edges.end()); +} + +template +bound create_bound_towards_minimum(edge_list& edges) { + if (edges.size() == 1) { + if (is_horizontal(edges.front())) { + reverse_horizontal(edges.front()); + } + bound bnd; + std::swap(bnd.edges, edges); + return bnd; + } + auto next_edge = edges.begin(); + auto edge = next_edge; + ++next_edge; + bool edge_is_horizontal = is_horizontal(*edge); + if (edge_is_horizontal) { + reverse_horizontal(*edge); + } + bool next_edge_is_horizontal; + bool y_increasing_before_last_horizontal = false; // assume false at start + + while (next_edge != edges.end()) { + next_edge_is_horizontal = is_horizontal(*next_edge); + if ((!next_edge_is_horizontal && !edge_is_horizontal && edge->bot == next_edge->bot)) { + break; + } + if (!next_edge_is_horizontal && edge_is_horizontal) { + if (y_increasing_before_last_horizontal && (next_edge->bot == edge->bot || next_edge->bot == edge->top)) { + break; + } + } else if (!y_increasing_before_last_horizontal && !edge_is_horizontal && next_edge_is_horizontal && + (edge->bot == next_edge->top || edge->bot == next_edge->bot)) { + y_increasing_before_last_horizontal = true; + } + edge_is_horizontal = next_edge_is_horizontal; + edge = next_edge; + if (edge_is_horizontal) { + reverse_horizontal(*edge); + } + ++next_edge; + } + bound bnd; + if (next_edge == edges.end()) { + std::swap(edges, bnd.edges); + } else { + bnd.edges.reserve(static_cast(std::distance(edges.begin(), next_edge))); + std::move(edges.begin(), next_edge, std::back_inserter(bnd.edges)); + edges.erase(edges.begin(), next_edge); + } + std::reverse(bnd.edges.begin(), bnd.edges.end()); + return bnd; +} + +template +bound create_bound_towards_maximum(edge_list& edges) { + if (edges.size() == 1) { + bound bnd; + std::swap(bnd.edges, edges); + return bnd; + } + auto next_edge = edges.begin(); + auto edge = next_edge; + ++next_edge; + bool edge_is_horizontal = is_horizontal(*edge); + bool next_edge_is_horizontal; + bool y_decreasing_before_last_horizontal = false; // assume false at start + + while (next_edge != edges.end()) { + next_edge_is_horizontal = is_horizontal(*next_edge); + if ((!next_edge_is_horizontal && !edge_is_horizontal && edge->top == next_edge->top)) { + break; + } + if (!next_edge_is_horizontal && edge_is_horizontal) { + if (y_decreasing_before_last_horizontal && (next_edge->top == edge->bot || next_edge->top == edge->top)) { + break; + } + } else if (!y_decreasing_before_last_horizontal && !edge_is_horizontal && next_edge_is_horizontal && + (edge->top == next_edge->top || edge->top == next_edge->bot)) { + y_decreasing_before_last_horizontal = true; + } + edge_is_horizontal = next_edge_is_horizontal; + edge = next_edge; + ++next_edge; + } + bound bnd; + if (next_edge == edges.end()) { + std::swap(bnd.edges, edges); + } else { + bnd.edges.reserve(static_cast(std::distance(edges.begin(), next_edge))); + std::move(edges.begin(), next_edge, std::back_inserter(bnd.edges)); + edges.erase(edges.begin(), next_edge); + } + return bnd; +} + +template +void fix_horizontals(bound& bnd) { + + auto edge_itr = bnd.edges.begin(); + auto next_itr = std::next(edge_itr); + if (next_itr == bnd.edges.end()) { + return; + } + if (is_horizontal(*edge_itr) && next_itr->bot != edge_itr->top) { + reverse_horizontal(*edge_itr); + } + auto prev_itr = edge_itr++; + while (edge_itr != bnd.edges.end()) { + if (is_horizontal(*edge_itr) && prev_itr->top != edge_itr->bot) { + reverse_horizontal(*edge_itr); + } + prev_itr = edge_itr; + ++edge_itr; + } +} + +template +void move_horizontals_on_left_to_right(bound& left_bound, bound& right_bound) { + // We want all the horizontal segments that are at the same Y as the minimum to be on the right + // bound + auto edge_itr = left_bound.edges.begin(); + while (edge_itr != left_bound.edges.end()) { + if (!is_horizontal(*edge_itr)) { + break; + } + reverse_horizontal(*edge_itr); + ++edge_itr; + } + if (edge_itr == left_bound.edges.begin()) { + return; + } + std::reverse(left_bound.edges.begin(), edge_itr); + auto dist = std::distance(left_bound.edges.begin(), edge_itr); + std::move(left_bound.edges.begin(), edge_itr, std::back_inserter(right_bound.edges)); + left_bound.edges.erase(left_bound.edges.begin(), edge_itr); + std::rotate(right_bound.edges.begin(), std::prev(right_bound.edges.end(), dist), right_bound.edges.end()); +} + +template +void add_ring_to_local_minima_list(edge_list& edges, local_minimum_list& minima_list, polygon_type poly_type) { + + if (edges.empty()) { + return; + } + // Adjust the order of the ring so we start on a local maximum + // therefore we start right away on a bound. + start_list_on_local_maximum(edges); + + bound_ptr first_minimum = nullptr; + bound_ptr last_maximum = nullptr; + while (!edges.empty()) { + interrupt_check(); // Check for interruptions + bool lm_minimum_has_horizontal = false; + auto to_minimum = create_bound_towards_minimum(edges); + if (edges.empty()) { + throw std::runtime_error("Edges is empty after only creating a single bound."); + } + auto to_maximum = create_bound_towards_maximum(edges); + fix_horizontals(to_minimum); + fix_horizontals(to_maximum); + auto to_max_first_non_horizontal = to_maximum.edges.begin(); + auto to_min_first_non_horizontal = to_minimum.edges.begin(); + bool minimum_is_left = true; + while (to_max_first_non_horizontal != to_maximum.edges.end() && is_horizontal(*to_max_first_non_horizontal)) { + lm_minimum_has_horizontal = true; + ++to_max_first_non_horizontal; + } + while (to_min_first_non_horizontal != to_minimum.edges.end() && is_horizontal(*to_min_first_non_horizontal)) { + lm_minimum_has_horizontal = true; + ++to_min_first_non_horizontal; + } + + if (to_max_first_non_horizontal == to_maximum.edges.end() || + to_min_first_non_horizontal == to_minimum.edges.end()) { + throw std::runtime_error("should not have a horizontal only bound for a ring"); + } + + if (lm_minimum_has_horizontal) { + if (to_max_first_non_horizontal->bot.x > to_min_first_non_horizontal->bot.x) { + minimum_is_left = true; + move_horizontals_on_left_to_right(to_minimum, to_maximum); + } else { + minimum_is_left = false; + move_horizontals_on_left_to_right(to_maximum, to_minimum); + } + } else { + if (to_max_first_non_horizontal->dx > to_min_first_non_horizontal->dx) { + minimum_is_left = false; + } else { + minimum_is_left = true; + } + } + assert(!to_minimum.edges.empty()); + assert(!to_maximum.edges.empty()); + auto const& min_front = to_minimum.edges.front(); + if (last_maximum) { + to_minimum.maximum_bound = last_maximum; + } + to_minimum.poly_type = poly_type; + to_maximum.poly_type = poly_type; + if (!minimum_is_left) { + to_minimum.side = edge_right; + to_maximum.side = edge_left; + to_minimum.winding_delta = -1; + to_maximum.winding_delta = 1; + minima_list.emplace_back(std::move(to_maximum), std::move(to_minimum), min_front.bot.y, + lm_minimum_has_horizontal); + if (!last_maximum) { + first_minimum = &(minima_list.back().right_bound); + } else { + last_maximum->maximum_bound = &(minima_list.back().right_bound); + } + last_maximum = &(minima_list.back().left_bound); + } else { + to_minimum.side = edge_left; + to_maximum.side = edge_right; + to_minimum.winding_delta = -1; + to_maximum.winding_delta = 1; + minima_list.emplace_back(std::move(to_minimum), std::move(to_maximum), min_front.bot.y, + lm_minimum_has_horizontal); + if (!last_maximum) { + first_minimum = &(minima_list.back().left_bound); + } else { + last_maximum->maximum_bound = &(minima_list.back().left_bound); + } + last_maximum = &(minima_list.back().right_bound); + } + } + last_maximum->maximum_bound = first_minimum; + first_minimum->maximum_bound = last_maximum; +} + +template +void initialize_lm(local_minimum_ptr_list_itr& lm) { + if (!(*lm)->left_bound.edges.empty()) { + (*lm)->left_bound.current_edge = (*lm)->left_bound.edges.begin(); + (*lm)->left_bound.next_edge = std::next((*lm)->left_bound.current_edge); + (*lm)->left_bound.current_x = static_cast((*lm)->left_bound.current_edge->bot.x); + (*lm)->left_bound.winding_count = 0; + (*lm)->left_bound.winding_count2 = 0; + (*lm)->left_bound.side = edge_left; + (*lm)->left_bound.ring = nullptr; + } + if (!(*lm)->right_bound.edges.empty()) { + (*lm)->right_bound.current_edge = (*lm)->right_bound.edges.begin(); + (*lm)->right_bound.next_edge = std::next((*lm)->right_bound.current_edge); + (*lm)->right_bound.current_x = static_cast((*lm)->right_bound.current_edge->bot.x); + (*lm)->right_bound.winding_count = 0; + (*lm)->right_bound.winding_count2 = 0; + (*lm)->right_bound.side = edge_right; + (*lm)->right_bound.ring = nullptr; + } +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/point.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/point.hpp new file mode 100644 index 000000000..1911914d8 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/point.hpp @@ -0,0 +1,110 @@ +#pragma once + +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct point; + +template +using point_ptr = point*; + +template +using const_point_ptr = point* const; + +template +struct ring; + +template +using ring_ptr = ring*; + +template +using const_ring_ptr = ring* const; + +template +struct point { + using coordinate_type = T; + ring_ptr ring; + T x; + T y; + point_ptr next; + point_ptr prev; + + point() : ring(nullptr), x(0), y(0), prev(this), next(this) { + } + + point(T x_, T y_) : ring(nullptr), x(x_), y(y_), next(this), prev(this) { + } + + point(ring_ptr ring_, mapbox::geometry::point const& pt) + : ring(ring_), x(pt.x), y(pt.y), next(this), prev(this) { + } + + point(ring_ptr ring_, mapbox::geometry::point const& pt, point_ptr before_this_point) + : ring(ring_), x(pt.x), y(pt.y), next(before_this_point), prev(before_this_point->prev) { + before_this_point->prev = this; + prev->next = this; + } +}; + +template +using point_vector = std::vector>; + +template +using point_vector_itr = typename point_vector::iterator; + +template +bool operator==(point const& lhs, point const& rhs) { + return lhs.x == rhs.x && lhs.y == rhs.y; +} + +template +bool operator==(mapbox::geometry::point const& lhs, point const& rhs) { + return lhs.x == rhs.x && lhs.y == rhs.y; +} + +template +bool operator==(point const& lhs, mapbox::geometry::point const& rhs) { + return lhs.x == rhs.x && lhs.y == rhs.y; +} + +template +bool operator!=(point const& lhs, point const& rhs) { + return lhs.x != rhs.x || lhs.y != rhs.y; +} + +template +bool operator!=(mapbox::geometry::point const& lhs, point const& rhs) { + return lhs.x != rhs.x || lhs.y != rhs.y; +} + +template +bool operator!=(point const& lhs, mapbox::geometry::point const& rhs) { + return lhs.x != rhs.x || lhs.y != rhs.y; +} + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, const point& p) { + out << " point at: " << p.x << ", " << p.y; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const mapbox::geometry::point& p) { + out << " point at: " << p.x << ", " << p.y; + return out; +} +#endif +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/process_horizontal.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/process_horizontal.hpp new file mode 100644 index 000000000..5dbbe54aa --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/process_horizontal.hpp @@ -0,0 +1,256 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +active_bound_list_itr process_horizontal_left_to_right(T scanline_y, + active_bound_list_itr& horz_bound, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + auto horizontal_itr_behind = horz_bound; + bool shifted = false; + bool is_maxima_edge = is_maxima(horz_bound, scanline_y); + auto bound_max_pair = active_bounds.end(); + if (is_maxima_edge) { + bound_max_pair = get_maxima_pair(horz_bound, active_bounds); + } + + auto hp_itr = rings.current_hp_itr; + while (hp_itr != rings.hot_pixels.end() && + (hp_itr->y > scanline_y || (hp_itr->y == scanline_y && hp_itr->x < (*horz_bound)->current_edge->bot.x))) { + ++hp_itr; + } + + auto bnd = std::next(horz_bound); + + while (bnd != active_bounds.end()) { + if (*bnd == nullptr) { + ++bnd; + continue; + } + // this code block inserts extra coords into horizontal edges (in output + // polygons) wherever hot pixels touch these horizontal edges. This helps + //'simplifying' polygons (ie if the Simplify property is set). + while (hp_itr != rings.hot_pixels.end() && hp_itr->y == scanline_y && + hp_itr->x < wround((*bnd)->current_x) && hp_itr->x < (*horz_bound)->current_edge->top.x) { + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + } + ++hp_itr; + } + + if ((*bnd)->current_x > static_cast((*horz_bound)->current_edge->top.x)) { + break; + } + + // Also break if we've got to the end of an intermediate horizontal edge ... + // nb: Smaller Dx's are to the right of larger Dx's ABOVE the horizontal. + if (wround((*bnd)->current_x) == (*horz_bound)->current_edge->top.x && + (*horz_bound)->next_edge != (*horz_bound)->edges.end() && + (*horz_bound)->current_edge->dx < (*horz_bound)->next_edge->dx) { + break; + } + + // note: may be done multiple times + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), mapbox::geometry::point(wround((*bnd)->current_x), scanline_y), + rings); + } + + // OK, so far we're still in range of the horizontal Edge but make sure + // we're at the last of consec. horizontals when matching with eMaxPair + if (is_maxima_edge && bnd == bound_max_pair) { + if ((*horz_bound)->ring) { + add_local_maximum_point(*(*horz_bound), *(*bound_max_pair), (*horz_bound)->current_edge->top, rings, + active_bounds); + } + *bound_max_pair = nullptr; + *horz_bound = nullptr; + if (!shifted) { + ++horizontal_itr_behind; + } + return horizontal_itr_behind; + } + + intersect_bounds(*(*horz_bound), *(*bnd), mapbox::geometry::point(wround((*bnd)->current_x), scanline_y), + cliptype, subject_fill_type, clip_fill_type, rings, active_bounds); + std::iter_swap(horz_bound, bnd); + horz_bound = bnd; + ++bnd; + shifted = true; + } // end while (bnd != active_bounds.end()) + + if ((*horz_bound)->ring) { + while (hp_itr != rings.hot_pixels.end() && hp_itr->y == scanline_y && + hp_itr->x < (*horz_bound)->current_edge->top.x) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + ++hp_itr; + } + } + + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), (*horz_bound)->current_edge->top, rings); + } + + if ((*horz_bound)->next_edge != (*horz_bound)->edges.end()) { + next_edge_in_bound(*(*horz_bound), scanbeam); + } else { + *horz_bound = nullptr; + } + if (!shifted) { + ++horizontal_itr_behind; + } + return horizontal_itr_behind; +} + +template +active_bound_list_itr process_horizontal_right_to_left(T scanline_y, + active_bound_list_itr& horz_bound_fwd, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + auto next_bnd_itr = std::next(horz_bound_fwd); + bool is_maxima_edge = is_maxima(horz_bound_fwd, scanline_y); + auto bound_max_pair = active_bounds.rend(); + if (is_maxima_edge) { + bound_max_pair = active_bound_list_rev_itr(get_maxima_pair(horz_bound_fwd, active_bounds)); + --bound_max_pair; + } + auto hp_itr_fwd = rings.current_hp_itr; + while (hp_itr_fwd != rings.hot_pixels.end() && + (hp_itr_fwd->y < scanline_y || + (hp_itr_fwd->y == scanline_y && hp_itr_fwd->x < (*horz_bound_fwd)->current_edge->top.x))) { + ++hp_itr_fwd; + } + auto hp_itr = hot_pixel_rev_itr(hp_itr_fwd); + + auto bnd = active_bound_list_rev_itr(horz_bound_fwd); + auto horz_bound = std::prev(bnd); + while (bnd != active_bounds.rend()) { + if (*bnd == nullptr) { + ++bnd; + continue; + } + // this code block inserts extra coords into horizontal edges (in output + // polygons) wherever hot pixels touch these horizontal edges. + while (hp_itr != rings.hot_pixels.rend() && hp_itr->y == scanline_y && + hp_itr->x > wround((*bnd)->current_x) && hp_itr->x > (*horz_bound)->current_edge->top.x) { + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + } + ++hp_itr; + } + + if ((*bnd)->current_x < static_cast((*horz_bound)->current_edge->top.x)) { + break; + } + + // Also break if we've got to the end of an intermediate horizontal edge ... + // nb: Smaller Dx's are to the right of larger Dx's ABOVE the horizontal. + if (wround((*bnd)->current_x) == (*horz_bound)->current_edge->top.x && + (*horz_bound)->next_edge != (*horz_bound)->edges.end() && + (*horz_bound)->current_edge->dx < (*horz_bound)->next_edge->dx) { + break; + } + + // note: may be done multiple times + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), mapbox::geometry::point(wround((*bnd)->current_x), scanline_y), + rings); + } + + // OK, so far we're still in range of the horizontal Edge but make sure + // we're at the last of consec. horizontals when matching with eMaxPair + if (is_maxima_edge && bnd == bound_max_pair) { + if ((*horz_bound)->ring) { + add_local_maximum_point(*(*horz_bound), *(*bound_max_pair), (*horz_bound)->current_edge->top, rings, + active_bounds); + } + *bound_max_pair = nullptr; + *horz_bound = nullptr; + return next_bnd_itr; + } + + intersect_bounds(*(*bnd), *(*horz_bound), mapbox::geometry::point(wround((*bnd)->current_x), scanline_y), + cliptype, subject_fill_type, clip_fill_type, rings, active_bounds); + std::iter_swap(horz_bound, bnd); + horz_bound = bnd; + ++bnd; + } // end while (bnd != active_bounds.rend()) + + if ((*horz_bound)->ring) { + while (hp_itr != rings.hot_pixels.rend() && hp_itr->y == scanline_y && + hp_itr->x > (*horz_bound)->current_edge->top.x) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + ++hp_itr; + } + } + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), (*horz_bound)->current_edge->top, rings); + } + + if ((*horz_bound)->next_edge != (*horz_bound)->edges.end()) { + next_edge_in_bound(*(*horz_bound), scanbeam); + } else { + *horz_bound = nullptr; + } + return next_bnd_itr; +} + +template +active_bound_list_itr process_horizontal(T scanline_y, + active_bound_list_itr& horz_bound, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + if ((*horz_bound)->current_edge->bot.x < (*horz_bound)->current_edge->top.x) { + return process_horizontal_left_to_right(scanline_y, horz_bound, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + } else { + return process_horizontal_right_to_left(scanline_y, horz_bound, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + } +} + +template +void process_horizontals(T scanline_y, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + for (auto bnd_itr = active_bounds.begin(); bnd_itr != active_bounds.end();) { + if (*bnd_itr != nullptr && current_edge_is_horizontal(bnd_itr)) { + bnd_itr = process_horizontal(scanline_y, bnd_itr, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + } else { + ++bnd_itr; + } + } + active_bounds.erase(std::remove(active_bounds.begin(), active_bounds.end(), nullptr), active_bounds.end()); +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/process_maxima.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/process_maxima.hpp new file mode 100644 index 000000000..17a9aaaaa --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/process_maxima.hpp @@ -0,0 +1,123 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +active_bound_list_itr do_maxima(active_bound_list_itr& bnd, + active_bound_list_itr& bndMaxPair, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& manager, + active_bound_list& active_bounds) { + auto bnd_next = std::next(bnd); + auto return_bnd = bnd; + bool skipped = false; + while (bnd_next != active_bounds.end() && bnd_next != bndMaxPair) { + if (*bnd_next == nullptr) { + ++bnd_next; + continue; + } + skipped = true; + intersect_bounds(*(*bnd), *(*bnd_next), (*bnd)->current_edge->top, cliptype, subject_fill_type, clip_fill_type, + manager, active_bounds); + std::iter_swap(bnd, bnd_next); + bnd = bnd_next; + ++bnd_next; + } + + if ((*bnd)->ring && (*bndMaxPair)->ring) { + add_local_maximum_point(*(*bnd), *(*bndMaxPair), (*bnd)->current_edge->top, manager, active_bounds); + } else if ((*bnd)->ring || (*bndMaxPair)->ring) { + throw std::runtime_error("DoMaxima error"); + } + *bndMaxPair = nullptr; + *bnd = nullptr; + if (!skipped) { + ++return_bnd; + } + return return_bnd; +} + +template +void process_edges_at_top_of_scanbeam(T top_y, + active_bound_list& active_bounds, + scanbeam_list& scanbeam, + local_minimum_ptr_list const& minima_sorted, + local_minimum_ptr_list_itr& current_lm, + ring_manager& manager, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + for (auto bnd = active_bounds.begin(); bnd != active_bounds.end();) { + interrupt_check(); // Check for interruptions + if (*bnd == nullptr) { + ++bnd; + continue; + } + // 1. Process maxima, treating them as if they are "bent" horizontal edges, + // but exclude maxima with horizontal edges. + + bool is_maxima_edge = is_maxima(bnd, top_y); + + if (is_maxima_edge) { + auto bnd_max_pair = get_maxima_pair(bnd, active_bounds); + is_maxima_edge = ((bnd_max_pair == active_bounds.end() || !current_edge_is_horizontal(bnd_max_pair)) && + is_maxima(bnd_max_pair, top_y)); + if (is_maxima_edge) { + bnd = do_maxima(bnd, bnd_max_pair, cliptype, subject_fill_type, clip_fill_type, manager, active_bounds); + continue; + } + } + + // 2. Promote horizontal edges. + if (is_intermediate(bnd, top_y) && next_edge_is_horizontal(bnd)) { + if ((*bnd)->ring) { + insert_hot_pixels_in_path(*(*bnd), (*bnd)->current_edge->top, manager, false); + } + next_edge_in_bound(*(*bnd), scanbeam); + if ((*bnd)->ring) { + add_point_to_ring(*(*bnd), (*bnd)->current_edge->bot, manager); + } + } else { + (*bnd)->current_x = get_current_x(*((*bnd)->current_edge), top_y); + } + ++bnd; + } + active_bounds.erase(std::remove(active_bounds.begin(), active_bounds.end(), nullptr), active_bounds.end()); + + insert_horizontal_local_minima_into_ABL(top_y, minima_sorted, current_lm, active_bounds, manager, scanbeam, + cliptype, subject_fill_type, clip_fill_type); + + process_horizontals(top_y, active_bounds, manager, scanbeam, cliptype, subject_fill_type, clip_fill_type); + + // 4. Promote intermediate vertices + + for (auto bnd = active_bounds.begin(); bnd != active_bounds.end(); ++bnd) { + if (is_intermediate(bnd, top_y)) { + if ((*bnd)->ring) { + add_point_to_ring(*(*bnd), (*bnd)->current_edge->top, manager); + } + next_edge_in_bound(*(*bnd), scanbeam); + } + } +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/quick_clip.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/quick_clip.hpp new file mode 100644 index 000000000..1b137b976 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/quick_clip.hpp @@ -0,0 +1,139 @@ +#pragma once + +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { +namespace quick_clip { + +template +mapbox::geometry::point intersect(mapbox::geometry::point a, + mapbox::geometry::point b, + size_t edge, + mapbox::geometry::box const& box) { + switch (edge) { + case 0: + return mapbox::geometry::point( + mapbox::geometry::wagyu::wround(static_cast(a.x) + static_cast(b.x - a.x) * + static_cast(box.min.y - a.y) / + static_cast(b.y - a.y)), + box.min.y); + + case 1: + return mapbox::geometry::point( + box.max.x, + mapbox::geometry::wagyu::wround(static_cast(a.y) + static_cast(b.y - a.y) * + static_cast(box.max.x - a.x) / + static_cast(b.x - a.x))); + + case 2: + return mapbox::geometry::point( + mapbox::geometry::wagyu::wround(static_cast(a.x) + static_cast(b.x - a.x) * + static_cast(box.max.y - a.y) / + static_cast(b.y - a.y)), + box.max.y); + + default: // case 3 + return mapbox::geometry::point( + box.min.x, + mapbox::geometry::wagyu::wround(static_cast(a.y) + static_cast(b.y - a.y) * + static_cast(box.min.x - a.x) / + static_cast(b.x - a.x))); + } +} + +template +bool inside(mapbox::geometry::point p, size_t edge, mapbox::geometry::box const& b) { + switch (edge) { + case 0: + return p.y > b.min.y; + + case 1: + return p.x < b.max.x; + + case 2: + return p.y < b.max.y; + + default: // case 3 + return p.x > b.min.x; + } +} + +template +mapbox::geometry::linear_ring quick_lr_clip(mapbox::geometry::linear_ring const& ring, + mapbox::geometry::box const& b) { + mapbox::geometry::linear_ring out = ring; + + for (size_t edge = 0; edge < 4; edge++) { + if (out.size() > 0) { + mapbox::geometry::linear_ring in = out; + mapbox::geometry::point S = in[in.size() - 1]; + out.resize(0); + + for (size_t e = 0; e < in.size(); e++) { + mapbox::geometry::point E = in[e]; + + if (inside(E, edge, b)) { + if (!inside(S, edge, b)) { + out.push_back(intersect(S, E, edge, b)); + } + out.push_back(E); + } else if (inside(S, edge, b)) { + out.push_back(intersect(S, E, edge, b)); + } + + S = E; + } + } + } + + if (out.size() < 3) { + out.clear(); + return out; + } + // Close the ring if the first/last point was outside + if (out[0] != out[out.size() - 1]) { + out.push_back(out[0]); + } + return out; +} +} // namespace quick_clip + +template +mapbox::geometry::multi_polygon +clip(mapbox::geometry::polygon const& poly, mapbox::geometry::box const& b, fill_type subject_fill_type) { + mapbox::geometry::multi_polygon result; + wagyu clipper; + for (auto const& lr : poly) { + auto new_lr = quick_clip::quick_lr_clip(lr, b); + if (!new_lr.empty()) { + clipper.add_ring(new_lr, polygon_type_subject); + } + } + clipper.execute(clip_type_union, result, subject_fill_type, fill_type_even_odd); + return result; +} + +template +mapbox::geometry::multi_polygon +clip(mapbox::geometry::multi_polygon const& mp, mapbox::geometry::box const& b, fill_type subject_fill_type) { + mapbox::geometry::multi_polygon result; + wagyu clipper; + for (auto const& poly : mp) { + for (auto const& lr : poly) { + auto new_lr = quick_clip::quick_lr_clip(lr, b); + if (!new_lr.empty()) { + clipper.add_ring(new_lr, polygon_type_subject); + } + } + } + clipper.execute(clip_type_union, result, subject_fill_type, fill_type_even_odd); + return result; +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/ring.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/ring.hpp new file mode 100644 index 000000000..eb4ade7f0 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/ring.hpp @@ -0,0 +1,633 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEBUG +#include +#include +#include +#include +// +// void* callstack[128]; +// int i, frames = backtrace(callstack, 128); +// char** strs = backtrace_symbols(callstack, frames); +// for (i = 0; i < frames; ++i) { +// printf("%s\n", strs[i]); +// } +// free(strs); +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +double area_from_point(point_ptr op, std::size_t& size, mapbox::geometry::box& bbox) { + point_ptr startOp = op; + size = 0; + double a = 0.0; + T min_x = op->x; + T max_x = op->x; + T min_y = op->y; + T max_y = op->y; + do { + ++size; + if (op->x > max_x) { + max_x = op->x; + } else if (op->x < min_x) { + min_x = op->x; + } + if (op->y > max_y) { + max_y = op->y; + } else if (op->y < min_y) { + min_y = op->y; + } + a += static_cast(op->prev->x + op->x) * static_cast(op->prev->y - op->y); + op = op->next; + } while (op != startOp); + bbox.min.x = min_x; + bbox.max.x = max_x; + bbox.min.y = min_y; + bbox.max.y = max_y; + return a * 0.5; +} + +// NOTE: ring and ring_ptr are forward declared in wagyu/point.hpp + +template +using ring_vector = std::vector>; + +template +struct ring { + std::size_t ring_index; // To support unset 0 is undefined and indexes offset by 1 + + std::size_t size_; // number of points in the ring + double area_; // area of the ring + mapbox::geometry::box bbox; // bounding box of the ring + + ring_ptr parent; + ring_vector children; + + point_ptr points; + point_ptr bottom_point; + bool is_hole_; + bool corrected; + + ring(ring const&) = delete; + ring& operator=(ring const&) = delete; + + ring() + : ring_index(0), + size_(0), + area_(std::numeric_limits::quiet_NaN()), + bbox({ 0, 0 }, { 0, 0 }), + parent(nullptr), + children(), + points(nullptr), + bottom_point(nullptr), + is_hole_(false), + corrected(false) { + } + + void reset_stats() { + area_ = std::numeric_limits::quiet_NaN(); + is_hole_ = false; + bbox.min.x = 0; + bbox.min.y = 0; + bbox.max.x = 0; + bbox.max.y = 0; + size_ = 0; + } + + void recalculate_stats() { + if (points != nullptr) { + area_ = area_from_point(points, size_, bbox); + is_hole_ = !(area_ > 0.0); + } + } + + void set_stats(double a, std::size_t s, mapbox::geometry::box const& b) { + bbox = b; + area_ = a; + size_ = s; + is_hole_ = !(area_ > 0.0); + } + + double area() { + if (std::isnan(area_)) { + recalculate_stats(); + } + return area_; + } + + bool is_hole() { + if (std::isnan(area_)) { + recalculate_stats(); + } + return is_hole_; + } + + std::size_t size() { + if (std::isnan(area_)) { + recalculate_stats(); + } + return size_; + } +}; + +template +using hot_pixel_vector = std::vector>; + +template +using hot_pixel_itr = typename hot_pixel_vector::iterator; + +template +using hot_pixel_rev_itr = typename hot_pixel_vector::reverse_iterator; + +template +struct ring_manager { + + ring_vector children; + point_vector all_points; + hot_pixel_vector hot_pixels; + hot_pixel_itr current_hp_itr; + std::deque> points; + std::deque> rings; + std::vector> storage; + std::size_t index; + + ring_manager(ring_manager const&) = delete; + ring_manager& operator=(ring_manager const&) = delete; + + ring_manager() + : children(), + all_points(), + hot_pixels(), + current_hp_itr(hot_pixels.end()), + points(), + rings(), + storage(), + index(0) { + } +}; + +template +void preallocate_point_memory(ring_manager& rings, std::size_t size) { + rings.storage.reserve(size); + rings.all_points.reserve(size); +} + +template +ring_ptr create_new_ring(ring_manager& manager) { + manager.rings.emplace_back(); + ring_ptr result = &manager.rings.back(); + result->ring_index = manager.index++; + return result; +} + +template +point_ptr create_new_point(ring_ptr r, mapbox::geometry::point const& pt, ring_manager& rings) { + point_ptr point; + if (rings.storage.size() < rings.storage.capacity()) { + rings.storage.emplace_back(r, pt); + point = &rings.storage.back(); + } else { + rings.points.emplace_back(r, pt); + point = &rings.points.back(); + } + rings.all_points.push_back(point); + return point; +} + +template +point_ptr create_new_point(ring_ptr r, + mapbox::geometry::point const& pt, + point_ptr before_this_point, + ring_manager& rings) { + point_ptr point; + if (rings.storage.size() < rings.storage.capacity()) { + rings.storage.emplace_back(r, pt, before_this_point); + point = &rings.storage.back(); + } else { + rings.points.emplace_back(r, pt, before_this_point); + point = &rings.points.back(); + } + rings.all_points.push_back(point); + return point; +} + +template +void set_to_children(ring_ptr r, ring_vector& children) { + for (auto& c : children) { + if (c == nullptr) { + c = r; + return; + } + } + children.push_back(r); +} + +template +void remove_from_children(ring_ptr r, ring_vector& children) { + for (auto& c : children) { + if (c == r) { + c = nullptr; + return; + } + } +} + +template +void assign_as_child(ring_ptr new_ring, ring_ptr parent, ring_manager& manager) { + // Assigning as a child assumes that this is + // a brand new ring. Therefore it does + // not have any existing relationships + if ((parent == nullptr && new_ring->is_hole()) || (parent != nullptr && new_ring->is_hole() == parent->is_hole())) { + throw std::runtime_error("Trying to assign a child that is the same orientation as the parent"); + } + auto& children = parent == nullptr ? manager.children : parent->children; + set_to_children(new_ring, children); + new_ring->parent = parent; +} + +template +void reassign_as_child(ring_ptr ring, ring_ptr parent, ring_manager& manager) { + // Reassigning a ring assumes it already + // has an existing parent + if ((parent == nullptr && ring->is_hole()) || (parent != nullptr && ring->is_hole() == parent->is_hole())) { + throw std::runtime_error("Trying to re-assign a child that is the same orientation as the parent"); + } + + // Remove the old child relationship + auto& old_children = ring->parent == nullptr ? manager.children : ring->parent->children; + remove_from_children(ring, old_children); + + // Add new child relationship + auto& children = parent == nullptr ? manager.children : parent->children; + set_to_children(ring, children); + ring->parent = parent; +} + +template +void assign_as_sibling(ring_ptr new_ring, ring_ptr sibling, ring_manager& manager) { + // Assigning as a sibling assumes that this is + // a brand new ring. Therefore it does + // not have any existing relationships + if (new_ring->is_hole() != sibling->is_hole()) { + throw std::runtime_error("Trying to assign to be a sibling that is not the same orientation as the sibling"); + } + auto& children = sibling->parent == nullptr ? manager.children : sibling->parent->children; + set_to_children(new_ring, children); + new_ring->parent = sibling->parent; +} + +template +void reassign_as_sibling(ring_ptr ring, ring_ptr sibling, ring_manager& manager) { + if (ring->parent == sibling->parent) { + return; + } + // Assigning as a sibling assumes that this is + // a brand new ring. Therefore it does + // not have any existing relationships + if (ring->is_hole() != sibling->is_hole()) { + throw std::runtime_error("Trying to assign to be a sibling that is not the same orientation as the sibling"); + } + // Remove the old child relationship + auto& old_children = ring->parent == nullptr ? manager.children : ring->parent->children; + remove_from_children(ring, old_children); + // Add new relationship + auto& children = sibling->parent == nullptr ? manager.children : sibling->parent->children; + set_to_children(ring, children); + ring->parent = sibling->parent; +} + +template +void ring1_replaces_ring2(ring_ptr ring1, ring_ptr ring2, ring_manager& manager) { + assert(ring1 != ring2); + auto& ring1_children = ring1 == nullptr ? manager.children : ring1->children; + for (auto& c : ring2->children) { + if (c == nullptr) { + continue; + } + c->parent = ring1; + set_to_children(c, ring1_children); + c = nullptr; + } + // Remove the old child relationship + auto& old_children = ring2->parent == nullptr ? manager.children : ring2->parent->children; + remove_from_children(ring2, old_children); + ring2->points = nullptr; + ring2->reset_stats(); +} + +template +void remove_points(point_ptr pt) { + if (pt != nullptr) { + pt->prev->next = nullptr; + while (pt != nullptr) { + point_ptr tmp = pt; + pt = pt->next; + tmp->next = nullptr; + tmp->prev = nullptr; + tmp->ring = nullptr; + } + } +} + +template +void remove_ring_and_points(ring_ptr r, + ring_manager& manager, + bool remove_children = true, + bool remove_from_parent = true) { + // Removes a ring and any children that might be + // under that ring. + for (auto& c : r->children) { + if (c == nullptr) { + continue; + } + if (remove_children) { + remove_ring_and_points(c, manager, true, false); + } + c = nullptr; + } + if (remove_from_parent) { + // Remove the old child relationship + auto& old_children = r->parent == nullptr ? manager.children : r->parent->children; + remove_from_children(r, old_children); + } + point_ptr pt = r->points; + if (pt != nullptr) { + pt->prev->next = nullptr; + while (pt != nullptr) { + point_ptr tmp = pt; + pt = pt->next; + tmp->next = nullptr; + tmp->prev = nullptr; + tmp->ring = nullptr; + } + } + r->points = nullptr; + r->reset_stats(); +} + +template +void remove_ring(ring_ptr r, ring_manager& manager, bool remove_children = true, bool remove_from_parent = true) { + // Removes a ring and any children that might be + // under that ring. + for (auto& c : r->children) { + if (c == nullptr) { + continue; + } + if (remove_children) { + remove_ring(c, manager, true, false); + } + c = nullptr; + } + if (remove_from_parent) { + // Remove the old child relationship + auto& old_children = r->parent == nullptr ? manager.children : r->parent->children; + remove_from_children(r, old_children); + } + r->points = nullptr; + r->reset_stats(); +} + +template +inline std::size_t ring_depth(ring_ptr r) { + std::size_t depth = 0; + if (!r) { + return depth; + } + while (r->parent) { + depth++; + r = r->parent; + } + return depth; +} + +template +inline bool ring_is_hole(ring_ptr r) { + // This is different then the "normal" way of determing if + // a ring is a hole or not because it uses the depth of the + // the ring to determine if it is a hole or not. This is only done + // intially when rings are output from Vatti. + return ring_depth(r) & 1; +} + +template +void set_next(const_point_ptr& node, const const_point_ptr& next_node) { + node->next = next_node; +} + +template +point_ptr get_next(const_point_ptr& node) { + return node->next; +} + +template +point_ptr get_prev(const_point_ptr& node) { + return node->prev; +} + +template +void set_prev(const_point_ptr& node, const const_point_ptr& prev_node) { + node->prev = prev_node; +} + +template +void init(const_point_ptr& node) { + set_next(node, node); + set_prev(node, node); +} + +template +void link_before(point_ptr& node, point_ptr& new_node) { + point_ptr prev_node = get_prev(node); + set_prev(new_node, prev_node); + set_next(new_node, node); + set_prev(node, new_node); + set_next(prev_node, new_node); +} + +template +void link_after(point_ptr& node, point_ptr& new_node) { + point_ptr next_node = get_next(node); + set_prev(new_node, node); + set_next(new_node, next_node); + set_next(node, new_node); + set_prev(next_node, new_node); +} + +template +void transfer_point(point_ptr& p, point_ptr& b, point_ptr& e) { + if (b != e) { + point_ptr prev_p = get_prev(p); + point_ptr prev_b = get_prev(b); + point_ptr prev_e = get_prev(e); + set_next(prev_e, p); + set_prev(p, prev_e); + set_next(prev_b, e); + set_prev(e, prev_b); + set_next(prev_p, b); + set_prev(b, prev_p); + } else { + link_before(p, b); + } +} + +template +void reverse_ring(point_ptr pp) { + if (!pp) { + return; + } + point_ptr pp1; + point_ptr pp2; + pp1 = pp; + do { + pp2 = pp1->next; + pp1->next = pp1->prev; + pp1->prev = pp2; + pp1 = pp2; + } while (pp1 != pp); +} + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, ring& r) { + out << " ring_index: " << r.ring_index << std::endl; + if (!r.parent) { + // out << " parent_ring ptr: nullptr" << std::endl; + out << " parent_index: -----" << std::endl; + } else { + // out << " parent_ring ptr: " << r.parent << std::endl; + out << " parent_ring idx: " << r.parent->ring_index << std::endl; + } + ring_ptr n = const_cast>(&r); + if (ring_is_hole(n)) { + out << " is_hole: true " << std::endl; + } else { + out << " is_hole: false " << std::endl; + } + auto pt_itr = r.points; + if (pt_itr) { + out << " area: " << r.area() << std::endl; + out << " points:" << std::endl; + out << " [[[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + while (pt_itr != r.points) { + out << "[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + } + out << "[" << pt_itr->x << "," << pt_itr->y << "]]]" << std::endl; + } else { + out << " area: NONE" << std::endl; + out << " points: NONE" << std::endl; + } + return out; +} + +template +std::string debug_ring_addresses(ring_ptr r) { + std::ostringstream out; + out << "Ring: " << r->ring_index << std::endl; + if (r->points == nullptr) { + out << " Ring has no points" << std::endl; + return out.str(); + } + auto pt_itr = r->points; + do { + out << " [" << pt_itr->x << "," << pt_itr->y << "] - " << pt_itr << std::endl; + pt_itr = pt_itr->next; + } while (pt_itr != r->points); + return out.str(); +} + +template +std::string output_as_polygon(ring_ptr r) { + std::ostringstream out; + + auto pt_itr = r->points; + if (pt_itr) { + out << "["; + out << "[[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + while (pt_itr != r->points) { + out << "[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + } + out << "[" << pt_itr->x << "," << pt_itr->y << "]]"; + for (auto const& c : r->children) { + if (c == nullptr) { + continue; + } + pt_itr = c->points; + if (pt_itr) { + out << ",[[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + while (pt_itr != c->points) { + out << "[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + } + out << "[" << pt_itr->x << "," << pt_itr->y << "]]"; + } + } + out << "]" << std::endl; + } else { + out << "[]" << std::endl; + } + + return out.str(); +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, ring_vector& rings) { + out << "START RING VECTOR" << std::endl; + for (auto& r : rings) { + if (r == nullptr || !r->points) { + continue; + } + out << " ring: " << r->ring_index << " - " << r << std::endl; + out << *r; + } + out << "END RING VECTOR" << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + std::deque>& rings) { + out << "START RING VECTOR" << std::endl; + for (auto& r : rings) { + if (!r.points) { + continue; + } + out << " ring: " << r.ring_index << std::endl; + out << r; + } + out << "END RING VECTOR" << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + hot_pixel_vector& hp_vec) { + out << "Hot Pixels: " << std::endl; + for (auto& hp : hp_vec) { + out << hp << std::endl; + } + return out; +} +#endif +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/ring_util.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/ring_util.hpp new file mode 100644 index 000000000..8aa39509e --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/ring_util.hpp @@ -0,0 +1,832 @@ +#pragma once + +#ifdef DEBUG +#include +// Example debug print for backtrace - only works on IOS +#include +#include +// +// void* callstack[128]; +// int i, frames = backtrace(callstack, 128); +// char** strs = backtrace_symbols(callstack, frames); +// for (i = 0; i < frames; ++i) { +// printf("%s\n", strs[i]); +// } +// free(strs); +#endif + +#include + +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +void set_hole_state(bound& bnd, active_bound_list const& active_bounds, ring_manager& rings) { + auto bnd_itr = std::find(active_bounds.rbegin(), active_bounds.rend(), &bnd); + ++bnd_itr; + bound_ptr bndTmp = nullptr; + // Find first non line ring to the left of current bound. + while (bnd_itr != active_bounds.rend()) { + if (*bnd_itr == nullptr) { + ++bnd_itr; + continue; + } + if ((*bnd_itr)->ring) { + if (!bndTmp) { + bndTmp = (*bnd_itr); + } else if (bndTmp->ring == (*bnd_itr)->ring) { + bndTmp = nullptr; + } + } + ++bnd_itr; + } + if (!bndTmp) { + bnd.ring->parent = nullptr; + rings.children.push_back(bnd.ring); + } else { + bnd.ring->parent = bndTmp->ring; + bndTmp->ring->children.push_back(bnd.ring); + } +} + +template +void update_current_hp_itr(T scanline_y, ring_manager& rings) { + while (rings.current_hp_itr->y > scanline_y) { + ++rings.current_hp_itr; + } +} + +template +struct hot_pixel_sorter { + inline bool operator()(mapbox::geometry::point const& pt1, mapbox::geometry::point const& pt2) { + if (pt1.y == pt2.y) { + return pt1.x < pt2.x; + } else { + return pt1.y > pt2.y; + } + } +}; + +// Due to the nature of floating point calculations +// and the high likely hood of values around X.5, we +// need to fudge what is X.5 some for our rounding. +const double rounding_offset = 1e-12; +const double rounding_offset_y = 5e-13; + +template +T round_towards_min(double val) { + // 0.5 rounds to 0 + // 0.0 rounds to 0 + // -0.5 rounds to -1 + return static_cast(std::ceil(val - 0.5 + rounding_offset)); +} + +template +T round_towards_max(double val) { + // 0.5 rounds to 1 + // 0.0 rounds to 0 + // -0.5 rounds to 0 + return static_cast(std::floor(val + 0.5 + rounding_offset)); +} + +template +inline T get_edge_min_x(edge const& edge, const T current_y) { + if (is_horizontal(edge)) { + if (edge.bot.x < edge.top.x) { + return edge.bot.x; + } else { + return edge.top.x; + } + } else if (edge.dx > 0.0) { + if (current_y == edge.top.y) { + return edge.top.x; + } else { + double lower_range_y = static_cast(current_y - edge.bot.y) - 0.5; + double return_val = static_cast(edge.bot.x) + edge.dx * lower_range_y; + T value = round_towards_min(return_val); + return value; + } + } else { + if (current_y == edge.bot.y) { + return edge.bot.x; + } else { + double return_val = static_cast(edge.bot.x) + + edge.dx * (static_cast(current_y - edge.bot.y) + 0.5 - rounding_offset_y); + T value = round_towards_min(return_val); + return value; + } + } +} + +template +inline T get_edge_max_x(edge const& edge, const T current_y) { + if (is_horizontal(edge)) { + if (edge.bot.x > edge.top.x) { + return edge.bot.x; + } else { + return edge.top.x; + } + } else if (edge.dx < 0.0) { + if (current_y == edge.top.y) { + return edge.top.x; + } else { + double lower_range_y = static_cast(current_y - edge.bot.y) - 0.5; + double return_val = static_cast(edge.bot.x) + edge.dx * lower_range_y; + T value = round_towards_max(return_val); + return value; + } + } else { + if (current_y == edge.bot.y) { + return edge.bot.x; + } else { + double return_val = static_cast(edge.bot.x) + + edge.dx * (static_cast(current_y - edge.bot.y) + 0.5 - rounding_offset_y); + T value = round_towards_max(return_val); + return value; + } + } +} + +template +void hot_pixel_set_left_to_right(T y, + T start_x, + T end_x, + bound& bnd, + ring_manager& rings, + hot_pixel_itr& itr, + hot_pixel_itr& end, + bool add_end_point) { + T x_min = get_edge_min_x(*(bnd.current_edge), y); + x_min = std::max(x_min, start_x); + T x_max = get_edge_max_x(*(bnd.current_edge), y); + x_max = std::min(x_max, end_x); + for (; itr != end; ++itr) { + if (itr->x < x_min) { + continue; + } + if (itr->x > x_max) { + break; + } + if (!add_end_point && itr->x == end_x) { + continue; + } + point_ptr op = bnd.ring->points; + bool to_front = (bnd.side == edge_left); + if (to_front && (*itr == *op)) { + continue; + } else if (!to_front && (*itr == *op->prev)) { + continue; + } + point_ptr new_point = create_new_point(bnd.ring, *itr, op, rings); + if (to_front) { + bnd.ring->points = new_point; + } + } +} + +template +void hot_pixel_set_right_to_left(T y, + T start_x, + T end_x, + bound& bnd, + ring_manager& rings, + hot_pixel_rev_itr& itr, + hot_pixel_rev_itr& end, + bool add_end_point) { + T x_min = get_edge_min_x(*(bnd.current_edge), y); + x_min = std::max(x_min, end_x); + T x_max = get_edge_max_x(*(bnd.current_edge), y); + x_max = std::min(x_max, start_x); + for (; itr != end; ++itr) { + if (itr->x > x_max) { + continue; + } + if (itr->x < x_min) { + break; + } + if (!add_end_point && itr->x == end_x) { + continue; + } + point_ptr op = bnd.ring->points; + bool to_front = (bnd.side == edge_left); + if (to_front && (*itr == *op)) { + continue; + } else if (!to_front && (*itr == *op->prev)) { + continue; + } + point_ptr new_point = create_new_point(bnd.ring, *itr, op, rings); + if (to_front) { + bnd.ring->points = new_point; + } + } +} + +template +void sort_hot_pixels(ring_manager& rings) { + std::sort(rings.hot_pixels.begin(), rings.hot_pixels.end(), hot_pixel_sorter()); + auto last = std::unique(rings.hot_pixels.begin(), rings.hot_pixels.end()); + rings.hot_pixels.erase(last, rings.hot_pixels.end()); +} + +template +void insert_hot_pixels_in_path(bound& bnd, + mapbox::geometry::point const& end_pt, + ring_manager& rings, + bool add_end_point) { + if (end_pt == bnd.last_point) { + return; + } + + T start_y = bnd.last_point.y; + T start_x = bnd.last_point.x; + T end_y = end_pt.y; + T end_x = end_pt.x; + + auto itr = rings.current_hp_itr; + while (itr->y <= start_y && itr != rings.hot_pixels.begin()) { + --itr; + } + if (start_x > end_x) { + for (; itr != rings.hot_pixels.end();) { + if (itr->y > start_y) { + ++itr; + continue; + } + if (itr->y < end_y) { + break; + } + T y = itr->y; + auto last_itr = hot_pixel_rev_itr(itr); + while (itr != rings.hot_pixels.end() && itr->y == y) { + ++itr; + } + auto first_itr = hot_pixel_rev_itr(itr); + bool add_end_point_itr = (y != end_pt.y || add_end_point); + hot_pixel_set_right_to_left(y, start_x, end_x, bnd, rings, first_itr, last_itr, add_end_point_itr); + } + } else { + for (; itr != rings.hot_pixels.end();) { + if (itr->y > start_y) { + ++itr; + continue; + } + if (itr->y < end_y) { + break; + } + T y = itr->y; + auto first_itr = itr; + while (itr != rings.hot_pixels.end() && itr->y == y) { + ++itr; + } + auto last_itr = itr; + bool add_end_point_itr = (y != end_pt.y || add_end_point); + hot_pixel_set_left_to_right(y, start_x, end_x, bnd, rings, first_itr, last_itr, add_end_point_itr); + } + } + bnd.last_point = end_pt; +} + +template +void add_to_hot_pixels(mapbox::geometry::point const& pt, ring_manager& rings) { + rings.hot_pixels.push_back(pt); +} + +template +void add_first_point(bound& bnd, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + + ring_ptr r = create_new_ring(rings); + bnd.ring = r; + r->points = create_new_point(r, pt, rings); + set_hole_state(bnd, active_bounds, rings); + bnd.last_point = pt; +} + +template +void add_point_to_ring(bound& bnd, mapbox::geometry::point const& pt, ring_manager& rings) { + assert(bnd.ring); + // Handle hot pixels + insert_hot_pixels_in_path(bnd, pt, rings, false); + + // bnd.ring->points is the 'Left-most' point & bnd.ring->points->prev is the + // 'Right-most' + point_ptr op = bnd.ring->points; + bool to_front = (bnd.side == edge_left); + if (to_front && (pt == *op)) { + return; + } else if (!to_front && (pt == *op->prev)) { + return; + } + point_ptr new_point = create_new_point(bnd.ring, pt, bnd.ring->points, rings); + if (to_front) { + bnd.ring->points = new_point; + } +} + +template +void add_point(bound& bnd, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + if (bnd.ring == nullptr) { + add_first_point(bnd, active_bounds, pt, rings); + } else { + add_point_to_ring(bnd, pt, rings); + } +} + +template +void add_local_minimum_point(bound& b1, + bound& b2, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + if (is_horizontal(*b2.current_edge) || (b1.current_edge->dx > b2.current_edge->dx)) { + add_point(b1, active_bounds, pt, rings); + b2.last_point = pt; + b2.ring = b1.ring; + b1.side = edge_left; + b2.side = edge_right; + } else { + add_point(b2, active_bounds, pt, rings); + b1.last_point = pt; + b1.ring = b2.ring; + b1.side = edge_right; + b2.side = edge_left; + } +} + +template +inline double get_dx(point const& pt1, point const& pt2) { + if (pt1.y == pt2.y) { + return std::numeric_limits::infinity(); + } else { + return static_cast(pt2.x - pt1.x) / static_cast(pt2.y - pt1.y); + } +} + +template +bool first_is_bottom_point(const_point_ptr btmPt1, const_point_ptr btmPt2) { + point_ptr p = btmPt1->prev; + while ((*p == *btmPt1) && (p != btmPt1)) { + p = p->prev; + } + double dx1p = std::fabs(get_dx(*btmPt1, *p)); + + p = btmPt1->next; + while ((*p == *btmPt1) && (p != btmPt1)) { + p = p->next; + } + double dx1n = std::fabs(get_dx(*btmPt1, *p)); + + p = btmPt2->prev; + while ((*p == *btmPt2) && (p != btmPt2)) { + p = p->prev; + } + double dx2p = std::fabs(get_dx(*btmPt2, *p)); + + p = btmPt2->next; + while ((*p == *btmPt2) && (p != btmPt2)) { + p = p->next; + } + double dx2n = std::fabs(get_dx(*btmPt2, *p)); + + if (values_are_equal(std::max(dx1p, dx1n), std::max(dx2p, dx2n)) && + values_are_equal(std::min(dx1p, dx1n), std::min(dx2p, dx2n))) { + std::size_t s = 0; + mapbox::geometry::box bbox({ 0, 0 }, { 0, 0 }); + return area_from_point(btmPt1, s, bbox) > 0.0; // if otherwise identical use orientation + } else { + return (greater_than_or_equal(dx1p, dx2p) && greater_than_or_equal(dx1p, dx2n)) || + (greater_than_or_equal(dx1n, dx2p) && greater_than_or_equal(dx1n, dx2n)); + } +} + +template +point_ptr get_bottom_point(point_ptr pp) { + point_ptr dups = nullptr; + point_ptr p = pp->next; + while (p != pp) { + if (p->y > pp->y) { + pp = p; + dups = nullptr; + } else if (p->y == pp->y && p->x <= pp->x) { + if (p->x < pp->x) { + dups = nullptr; + pp = p; + } else { + if (p->next != pp && p->prev != pp) { + dups = p; + } + } + } + p = p->next; + } + if (dups) { + // there appears to be at least 2 vertices at bottom_point so ... + while (dups != p) { + if (!first_is_bottom_point(p, dups)) { + pp = dups; + } + dups = dups->next; + while (*dups != *pp) { + dups = dups->next; + } + } + } + return pp; +} + +template +ring_ptr get_lower_most_ring(ring_ptr ring1, ring_ptr ring2) { + // work out which polygon fragment has the correct hole state ... + if (!ring1->bottom_point) { + ring1->bottom_point = get_bottom_point(ring1->points); + } + if (!ring2->bottom_point) { + ring2->bottom_point = get_bottom_point(ring2->points); + } + point_ptr pt1 = ring1->bottom_point; + point_ptr pt2 = ring2->bottom_point; + if (pt1->y > pt2->y) { + return ring1; + } else if (pt1->y < pt2->y) { + return ring2; + } else if (pt1->x < pt2->x) { + return ring1; + } else if (pt1->x > pt2->x) { + return ring2; + } else if (pt1->next == pt1) { + return ring2; + } else if (pt2->next == pt2) { + return ring1; + } else if (first_is_bottom_point(pt1, pt2)) { + return ring1; + } else { + return ring2; + } +} + +template +bool ring1_child_below_ring2(ring_ptr ring1, ring_ptr ring2) { + do { + ring1 = ring1->parent; + if (ring1 == ring2) { + return true; + } + } while (ring1); + return false; +} + +template +void update_points_ring(ring_ptr ring) { + point_ptr op = ring->points; + do { + op->ring = ring; + op = op->prev; + } while (op != ring->points); +} + +template +void append_ring(bound& b1, bound& b2, active_bound_list& active_bounds, ring_manager& manager) { + // get the start and ends of both output polygons ... + ring_ptr outRec1 = b1.ring; + ring_ptr outRec2 = b2.ring; + + ring_ptr keep_ring; + bound_ptr keep_bound; + ring_ptr remove_ring; + bound_ptr remove_bound; + if (ring1_child_below_ring2(outRec1, outRec2)) { + keep_ring = outRec2; + keep_bound = &b2; + remove_ring = outRec1; + remove_bound = &b1; + } else if (ring1_child_below_ring2(outRec2, outRec1)) { + keep_ring = outRec1; + keep_bound = &b1; + remove_ring = outRec2; + remove_bound = &b2; + } else if (outRec1 == get_lower_most_ring(outRec1, outRec2)) { + keep_ring = outRec1; + keep_bound = &b1; + remove_ring = outRec2; + remove_bound = &b2; + } else { + keep_ring = outRec2; + keep_bound = &b2; + remove_ring = outRec1; + remove_bound = &b1; + } + + // get the start and ends of both output polygons and + // join b2 poly onto b1 poly and delete pointers to b2 ... + + point_ptr p1_lft = keep_ring->points; + point_ptr p1_rt = p1_lft->prev; + point_ptr p2_lft = remove_ring->points; + point_ptr p2_rt = p2_lft->prev; + + // join b2 poly onto b1 poly and delete pointers to b2 ... + if (keep_bound->side == edge_left) { + if (remove_bound->side == edge_left) { + // z y x a b c + reverse_ring(p2_lft); + p2_lft->next = p1_lft; + p1_lft->prev = p2_lft; + p1_rt->next = p2_rt; + p2_rt->prev = p1_rt; + keep_ring->points = p2_rt; + } else { + // x y z a b c + p2_rt->next = p1_lft; + p1_lft->prev = p2_rt; + p2_lft->prev = p1_rt; + p1_rt->next = p2_lft; + keep_ring->points = p2_lft; + } + } else { + if (remove_bound->side == edge_right) { + // a b c z y x + reverse_ring(p2_lft); + p1_rt->next = p2_rt; + p2_rt->prev = p1_rt; + p2_lft->next = p1_lft; + p1_lft->prev = p2_lft; + } else { + // a b c x y z + p1_rt->next = p2_lft; + p2_lft->prev = p1_rt; + p1_lft->prev = p2_rt; + p2_rt->next = p1_lft; + } + } + + keep_ring->bottom_point = nullptr; + bool keep_is_hole = ring_is_hole(keep_ring); + bool remove_is_hole = ring_is_hole(remove_ring); + + remove_ring->points = nullptr; + remove_ring->bottom_point = nullptr; + if (keep_is_hole != remove_is_hole) { + ring1_replaces_ring2(keep_ring->parent, remove_ring, manager); + } else { + ring1_replaces_ring2(keep_ring, remove_ring, manager); + } + + update_points_ring(keep_ring); + + // nb: safe because we only get here via AddLocalMaxPoly + keep_bound->ring = nullptr; + remove_bound->ring = nullptr; + + for (auto& b : active_bounds) { + if (b == nullptr) { + continue; + } + if (b->ring == remove_ring) { + b->ring = keep_ring; + b->side = keep_bound->side; + break; // Not sure why there is a break here but was transfered logic from angus + } + } +} + +template +void add_local_maximum_point(bound& b1, + bound& b2, + mapbox::geometry::point const& pt, + ring_manager& rings, + active_bound_list& active_bounds) { + insert_hot_pixels_in_path(b2, pt, rings, false); + add_point(b1, active_bounds, pt, rings); + if (b1.ring == b2.ring) { + b1.ring = nullptr; + b2.ring = nullptr; + // I am not certain that order is important here? + } else if (b1.ring->ring_index < b2.ring->ring_index) { + append_ring(b1, b2, active_bounds, rings); + } else { + append_ring(b2, b1, active_bounds, rings); + } +} + +enum point_in_polygon_result : std::int8_t { + point_on_polygon = -1, + point_inside_polygon = 0, + point_outside_polygon = 1 +}; + +template +point_in_polygon_result point_in_polygon(point const& pt, point_ptr op) { + // returns 0 if false, +1 if true, -1 if pt ON polygon boundary + point_in_polygon_result result = point_outside_polygon; + point_ptr startOp = op; + do { + if (op->next->y == pt.y) { + if ((op->next->x == pt.x) || (op->y == pt.y && ((op->next->x > pt.x) == (op->x < pt.x)))) { + return point_on_polygon; + } + } + if ((op->y < pt.y) != (op->next->y < pt.y)) { + if (op->x >= pt.x) { + if (op->next->x > pt.x) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } else { + double d = static_cast(op->x - pt.x) * static_cast(op->next->y - pt.y) - + static_cast(op->next->x - pt.x) * static_cast(op->y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0) == (op->next->y > op->y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } else { + if (op->next->x > pt.x) { + double d = static_cast(op->x - pt.x) * static_cast(op->next->y - pt.y) - + static_cast(op->next->x - pt.x) * static_cast(op->y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0) == (op->next->y > op->y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } + } + op = op->next; + } while (startOp != op); + return result; +} + +template +point_in_polygon_result point_in_polygon(mapbox::geometry::point const& pt, point_ptr op) { + // returns 0 if false, +1 if true, -1 if pt ON polygon boundary + point_in_polygon_result result = point_outside_polygon; + point_ptr startOp = op; + do { + double op_x = static_cast(op->x); + double op_y = static_cast(op->y); + double op_next_x = static_cast(op->next->x); + double op_next_y = static_cast(op->next->y); + if (values_are_equal(op_next_y, pt.y)) { + if (values_are_equal(op_next_x, pt.x) || + (values_are_equal(op_y, pt.y) && ((op_next_x > pt.x) == (op_x < pt.x)))) { + return point_on_polygon; + } + } + if ((op_y < pt.y) != (op_next_y < pt.y)) { + if (greater_than_or_equal(op_x, pt.x)) { + if (op_next_x > pt.x) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } else { + double d = (op_x - pt.x) * (op_next_y - pt.y) - (op_next_x - pt.x) * (op_y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0.0) == (op_next_y > op_y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } else { + if (op_next_x > pt.x) { + double d = (op_x - pt.x) * (op_next_y - pt.y) - (op_next_x - pt.x) * (op_y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0.0) == (op_next_y > op_y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } + } + op = op->next; + } while (startOp != op); + return result; +} + +template +bool is_convex(point_ptr edge) { + point_ptr prev = edge->prev; + point_ptr next = edge->next; + T v1x = edge->x - prev->x; + T v1y = edge->y - prev->y; + T v2x = next->x - edge->x; + T v2y = next->y - edge->y; + T cross = v1x * v2y - v2x * v1y; + if (cross < 0 && edge->ring->area() > 0) { + return true; + } else if (cross > 0 && edge->ring->area() < 0) { + return true; + } else { + return false; + } +} + +template +mapbox::geometry::point centroid_of_points(point_ptr edge) { + point_ptr prev = edge->prev; + point_ptr next = edge->next; + return { static_cast(prev->x + edge->x + next->x) / 3.0, + static_cast(prev->y + edge->y + next->y) / 3.0 }; +} + +template +point_in_polygon_result inside_or_outside_special(point_ptr first_pt, point_ptr other_poly) { + + // We are going to loop through all the points + // of the original triangle. The goal is to find a convex edge + // that with its next and previous forms a triangle with its centroid + // that is within the first ring. Then we will check the other polygon + // to see if it is within this polygon. + point_ptr itr = first_pt; + do { + if (is_convex(itr)) { + auto pt = centroid_of_points(itr); + if (point_inside_polygon == point_in_polygon(pt, first_pt)) { + return point_in_polygon(pt, other_poly); + } + } + itr = itr->next; + } while (itr != first_pt); + + throw std::runtime_error("Could not find a point within the polygon to test"); +} + +template +bool box2_contains_box1(mapbox::geometry::box const& box1, mapbox::geometry::box const& box2) { + return (box2.max.x >= box1.max.x && box2.max.y >= box1.max.y && box2.min.x <= box1.min.x && + box2.min.y <= box1.min.y); +} + +template +bool poly2_contains_poly1(ring_ptr ring1, ring_ptr ring2) { + if (!box2_contains_box1(ring1->bbox, ring2->bbox)) { + return false; + } + if (std::fabs(ring2->area()) < std::fabs(ring1->area())) { + return false; + } + point_ptr outpt1 = ring1->points->next; + point_ptr outpt2 = ring2->points->next; + point_ptr op = outpt1; + do { + // nb: PointInPolygon returns 0 if false, +1 if true, -1 if pt on polygon + point_in_polygon_result res = point_in_polygon(*op, outpt2); + if (res != point_on_polygon) { + return res == point_inside_polygon; + } + op = op->next; + } while (op != outpt1); + point_in_polygon_result res = inside_or_outside_special(outpt1, outpt2); + return res == point_inside_polygon; +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/scanbeam.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/scanbeam.hpp new file mode 100644 index 000000000..5f010abd1 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/scanbeam.hpp @@ -0,0 +1,45 @@ +#pragma once + +#include +#include + +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +using scanbeam_list = std::vector; + +template +void insert_sorted_scanbeam(scanbeam_list& scanbeam, T& t) { + typename scanbeam_list::iterator i = std::lower_bound(scanbeam.begin(), scanbeam.end(), t); + if (i == scanbeam.end() || t < *i) { + scanbeam.insert(i, t); + } +} + +template +bool pop_from_scanbeam(T& Y, scanbeam_list& scanbeam) { + if (scanbeam.empty()) { + return false; + } + + Y = scanbeam.back(); + scanbeam.pop_back(); + return true; +} + +template +void setup_scanbeam(local_minimum_list& minima_list, scanbeam_list& scanbeam) { + + scanbeam.reserve(minima_list.size()); + for (auto lm = minima_list.begin(); lm != minima_list.end(); ++lm) { + scanbeam.push_back(lm->y); + } + std::sort(scanbeam.begin(), scanbeam.end()); +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/snap_rounding.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/snap_rounding.hpp new file mode 100644 index 000000000..ea5cc8085 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/snap_rounding.hpp @@ -0,0 +1,191 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct hp_intersection_swap { + + ring_manager& manager; + + hp_intersection_swap(ring_manager& m) : manager(m) { + } + + void operator()(bound_ptr const& b1, bound_ptr const& b2) { + mapbox::geometry::point pt; + if (!get_edge_intersection(*(b1->current_edge), *(b2->current_edge), pt)) { + // LCOV_EXCL_START + throw std::runtime_error("Trying to find intersection of lines that do not intersect"); + // LCOV_EXCL_END + } + add_to_hot_pixels(round_point(pt), manager); + } +}; + +template +void process_hot_pixel_intersections(T top_y, active_bound_list& active_bounds, ring_manager& manager) { + if (active_bounds.empty()) { + return; + } + update_current_x(active_bounds, top_y); + bubble_sort(active_bounds.begin(), active_bounds.end(), intersection_compare(), + hp_intersection_swap(manager)); +} + +template +bool horizontals_at_top_scanbeam(T top_y, + active_bound_list_itr& bnd_curr, + active_bound_list& active_bounds, + ring_manager& manager) { + bool shifted = false; + auto& current_edge = (*bnd_curr)->current_edge; + (*bnd_curr)->current_x = static_cast(current_edge->top.x); + if (current_edge->bot.x < current_edge->top.x) { + // left to right + auto bnd_next = std::next(bnd_curr); + while (bnd_next != active_bounds.end() && + (*bnd_next == nullptr || (*bnd_next)->current_x < (*bnd_curr)->current_x)) { + if (*bnd_next != nullptr && (*bnd_next)->current_edge->top.y != top_y && + (*bnd_next)->current_edge->bot.y != top_y) { + mapbox::geometry::point pt(wround((*bnd_next)->current_x), top_y); + add_to_hot_pixels(pt, manager); + } + std::iter_swap(bnd_curr, bnd_next); + ++bnd_curr; + ++bnd_next; + shifted = true; + } + } else { + // right to left + if (bnd_curr != active_bounds.begin()) { + auto bnd_prev = std::prev(bnd_curr); + while (bnd_curr != active_bounds.begin() && + (*bnd_prev == nullptr || (*bnd_prev)->current_x > (*bnd_curr)->current_x)) { + if (*bnd_prev != nullptr && (*bnd_prev)->current_edge->top.y != top_y && + (*bnd_prev)->current_edge->bot.y != top_y) { + mapbox::geometry::point pt(wround((*bnd_prev)->current_x), top_y); + add_to_hot_pixels(pt, manager); + } + std::iter_swap(bnd_curr, bnd_prev); + --bnd_curr; + if (bnd_curr != active_bounds.begin()) { + --bnd_prev; + } + } + } + } + return shifted; +} + +template +void process_hot_pixel_edges_at_top_of_scanbeam(T top_y, + scanbeam_list& scanbeam, + active_bound_list& active_bounds, + ring_manager& manager) { + for (auto bnd = active_bounds.begin(); bnd != active_bounds.end();) { + if (*bnd == nullptr) { + ++bnd; + continue; + } + bound& current_bound = *(*bnd); + auto bnd_curr = bnd; + bool shifted = false; + auto& current_edge = current_bound.current_edge; + while (current_edge != current_bound.edges.end() && current_edge->top.y == top_y) { + add_to_hot_pixels(current_edge->top, manager); + if (is_horizontal(*current_edge)) { + if (horizontals_at_top_scanbeam(top_y, bnd_curr, active_bounds, manager)) { + shifted = true; + } + } + next_edge_in_bound(current_bound, scanbeam); + } + if (current_edge == current_bound.edges.end()) { + *bnd_curr = nullptr; + } + if (!shifted) { + ++bnd; + } + } + active_bounds.erase(std::remove(active_bounds.begin(), active_bounds.end(), nullptr), active_bounds.end()); +} + +template +void insert_local_minima_into_ABL_hot_pixel(T top_y, + local_minimum_ptr_list& minima_sorted, + local_minimum_ptr_list_itr& lm, + active_bound_list& active_bounds, + ring_manager& manager, + scanbeam_list& scanbeam) { + while (lm != minima_sorted.end() && (*lm)->y == top_y) { + add_to_hot_pixels((*lm)->left_bound.edges.front().bot, manager); + auto& left_bound = (*lm)->left_bound; + auto& right_bound = (*lm)->right_bound; + left_bound.current_edge = left_bound.edges.begin(); + left_bound.next_edge = std::next(left_bound.current_edge); + left_bound.current_x = static_cast(left_bound.current_edge->bot.x); + right_bound.current_edge = right_bound.edges.begin(); + right_bound.next_edge = std::next(right_bound.current_edge); + right_bound.current_x = static_cast(right_bound.current_edge->bot.x); + auto lb_abl_itr = insert_bound_into_ABL(left_bound, right_bound, active_bounds); + if (!current_edge_is_horizontal(lb_abl_itr)) { + insert_sorted_scanbeam(scanbeam, (*lb_abl_itr)->current_edge->top.y); + } + auto rb_abl_itr = std::next(lb_abl_itr); + if (!current_edge_is_horizontal(rb_abl_itr)) { + insert_sorted_scanbeam(scanbeam, (*rb_abl_itr)->current_edge->top.y); + } + ++lm; + } +} + +template +void build_hot_pixels(local_minimum_list& minima_list, ring_manager& manager) { + active_bound_list active_bounds; + scanbeam_list scanbeam; + T scanline_y = std::numeric_limits::max(); + + local_minimum_ptr_list minima_sorted; + minima_sorted.reserve(minima_list.size()); + for (auto& lm : minima_list) { + minima_sorted.push_back(&lm); + } + std::stable_sort(minima_sorted.begin(), minima_sorted.end(), local_minimum_sorter()); + local_minimum_ptr_list_itr current_lm = minima_sorted.begin(); + + setup_scanbeam(minima_list, scanbeam); + + // Estimate size for reserving hot pixels + std::size_t reserve = 0; + for (auto& lm : minima_list) { + reserve += lm.left_bound.edges.size() + 2; + reserve += lm.right_bound.edges.size() + 2; + } + manager.hot_pixels.reserve(reserve); + + while (pop_from_scanbeam(scanline_y, scanbeam) || current_lm != minima_sorted.end()) { + + process_hot_pixel_intersections(scanline_y, active_bounds, manager); + + insert_local_minima_into_ABL_hot_pixel(scanline_y, minima_sorted, current_lm, active_bounds, manager, scanbeam); + + process_hot_pixel_edges_at_top_of_scanbeam(scanline_y, scanbeam, active_bounds, manager); + } + preallocate_point_memory(manager, manager.hot_pixels.size()); + sort_hot_pixels(manager); +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/topology_correction.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/topology_correction.hpp new file mode 100644 index 000000000..4d7577829 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/topology_correction.hpp @@ -0,0 +1,1347 @@ +#pragma once + +#define _USE_MATH_DEFINES +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#ifdef DEBUG +#include +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct point_ptr_pair { + point_ptr op1; + point_ptr op2; + + constexpr point_ptr_pair(point_ptr o1, point_ptr o2) : op1(o1), op2(o2) { + } + + point_ptr_pair(point_ptr_pair const& p) = default; + + point_ptr_pair(point_ptr_pair&& p) : op1(std::move(p.op1)), op2(std::move(p.op2)) { + } + + point_ptr_pair& operator=(point_ptr_pair&& p) { + op1 = std::move(p.op1); + op2 = std::move(p.op2); + return *this; + } +}; + +#ifdef DEBUG + +template +inline std::basic_ostream& +operator<<(std::basic_ostream& out, + const std::unordered_multimap, point_ptr_pair>& dupe_ring) { + + out << " BEGIN CONNECTIONS: " << std::endl; + for (auto& r : dupe_ring) { + out << " Ring: "; + if (r.second.op1->ring) { + out << r.second.op1->ring->ring_index; + } else { + out << "---"; + } + out << " to "; + if (r.second.op2->ring) { + out << r.second.op2->ring->ring_index; + } else { + out << "---"; + } + out << " ( at " << r.second.op1->x << ", " << r.second.op1->y << " )"; + out << " Ring1 ( "; + if (r.second.op1->ring) { + out << "area: " << r.second.op1->ring->area << " parent: "; + if (r.second.op1->ring->parent) { + out << r.second.op1->ring->parent->ring_index; + } else { + out << "---"; + } + } else { + out << "---"; + } + out << " )"; + out << " Ring2 ( "; + if (r.second.op2->ring) { + out << "area: " << r.second.op2->ring->area << " parent: "; + if (r.second.op2->ring->parent) { + out << r.second.op2->ring->parent->ring_index; + } else { + out << "---"; + } + } else { + out << "---"; + } + out << " )"; + out << std::endl; + } + out << " END CONNECTIONS: " << std::endl; + return out; +} + +#endif + +template +bool find_intersect_loop(std::unordered_multimap, point_ptr_pair>& dupe_ring, + std::list, point_ptr_pair>>& iList, + ring_ptr ring_parent, + ring_ptr ring_origin, + ring_ptr ring_search, + std::set>& visited, + point_ptr orig_pt, + point_ptr prev_pt, + ring_manager& rings) { + { + auto range = dupe_ring.equal_range(ring_search); + // Check for direct connection + for (auto& it = range.first; it != range.second;) { + ring_ptr it_ring1 = it->second.op1->ring; + ring_ptr it_ring2 = it->second.op2->ring; + if (!it_ring1 || !it_ring2 || it_ring1 != ring_search || (!it_ring1->is_hole() && !it_ring2->is_hole())) { + it = dupe_ring.erase(it); + continue; + } + if (it_ring2 == ring_origin && (ring_parent == it_ring2 || ring_parent == it_ring2->parent) && + *prev_pt != *it->second.op2 && *orig_pt != *it->second.op2) { + iList.emplace_front(ring_search, it->second); + return true; + } + ++it; + } + } + auto range = dupe_ring.equal_range(ring_search); + visited.insert(ring_search); + // Check for connection through chain of other intersections + for (auto& it = range.first; it != range.second && it != dupe_ring.end() && it->first == ring_search; ++it) { + ring_ptr it_ring = it->second.op2->ring; + if (visited.count(it_ring) > 0 || it_ring == nullptr || + (ring_parent != it_ring && ring_parent != it_ring->parent) || value_is_zero(it_ring->area()) || + *prev_pt == *it->second.op2) { + continue; + } + if (find_intersect_loop(dupe_ring, iList, ring_parent, ring_origin, it_ring, visited, orig_pt, it->second.op2, + rings)) { + iList.emplace_front(ring_search, it->second); + return true; + } + } + return false; +} + +template +struct point_ptr_cmp { + inline bool operator()(point_ptr op1, point_ptr op2) { + if (op1->y != op2->y) { + return (op1->y > op2->y); + } else if (op1->x != op2->x) { + return (op1->x < op2->x); + } else { + std::size_t depth_1 = ring_depth(op1->ring); + std::size_t depth_2 = ring_depth(op2->ring); + return depth_1 > depth_2; + } + } +}; + +template +void correct_orientations(ring_manager& manager) { + for (auto& r : manager.rings) { + if (!r.points) { + continue; + } + r.recalculate_stats(); + if (r.size() < 3) { + remove_ring_and_points(&r, manager, false); + continue; + } + if (ring_is_hole(&r) != r.is_hole()) { + reverse_ring(r.points); + r.recalculate_stats(); + } + } +} + +template +point_vector sort_ring_points(ring_ptr r) { + point_vector sorted_points; + point_ptr point_itr = r->points; + point_ptr last_point = point_itr->prev; + while (point_itr != last_point) { + sorted_points.push_back(point_itr); + point_itr = point_itr->next; + } + sorted_points.push_back(last_point); + std::stable_sort(sorted_points.begin(), sorted_points.end(), [](point_ptr const& pt1, point_ptr const& pt2) { + if (pt1->y != pt2->y) { + return (pt1->y > pt2->y); + } + return (pt1->x < pt2->x); + }); + return sorted_points; +} + +template +ring_ptr correct_self_intersection(point_ptr pt1, point_ptr pt2, ring_manager& manager) { + if (pt1->ring != pt2->ring) { + return static_cast>(nullptr); + } + + ring_ptr ring = pt1->ring; + + // split the polygon into two ... + point_ptr pt3 = pt1->prev; + point_ptr pt4 = pt2->prev; + pt1->prev = pt4; + pt4->next = pt1; + pt2->prev = pt3; + pt3->next = pt2; + + ring_ptr new_ring = create_new_ring(manager); + std::size_t size_1 = 0; + std::size_t size_2 = 0; + mapbox::geometry::box box1({ 0, 0 }, { 0, 0 }); + mapbox::geometry::box box2({ 0, 0 }, { 0, 0 }); + double area_1 = area_from_point(pt1, size_1, box1); + double area_2 = area_from_point(pt2, size_2, box2); + + if (std::fabs(area_1) > std::fabs(area_2)) { + ring->points = pt1; + ring->set_stats(area_1, size_1, box1); + new_ring->points = pt2; + new_ring->set_stats(area_2, size_2, box2); + } else { + ring->points = pt2; + ring->set_stats(area_2, size_2, box2); + new_ring->points = pt1; + new_ring->set_stats(area_1, size_1, box1); + } + update_points_ring(new_ring); + return new_ring; +} + +template +void correct_repeated_points(ring_manager& manager, + ring_vector& new_rings, + point_vector_itr const& begin, + point_vector_itr const& end) { + for (auto itr1 = begin; itr1 != end; ++itr1) { + if ((*itr1)->ring == nullptr) { + continue; + } + for (auto itr2 = std::next(itr1); itr2 != end; ++itr2) { + if ((*itr2)->ring == nullptr) { + continue; + } + ring_ptr new_ring = correct_self_intersection(*itr1, *itr2, manager); + if (new_ring != nullptr) { + new_rings.push_back(new_ring); + } + } + } +} + +template +void find_and_correct_repeated_points(ring_ptr r, ring_manager& manager, ring_vector& new_rings) { + auto sorted_points = sort_ring_points(r); + // Find sets of repeated points + std::size_t count = 0; + auto prev_itr = sorted_points.begin(); + auto itr = std::next(prev_itr); + while (itr != sorted_points.end()) { + if (*(*prev_itr) == *(*(itr))) { + ++count; + ++prev_itr; + ++itr; + if (itr != sorted_points.end()) { + continue; + } else { + ++prev_itr; + } + } else { + ++prev_itr; + ++itr; + } + if (count == 0) { + continue; + } + auto first = prev_itr; + std::advance(first, -(static_cast(count) + 1)); + correct_repeated_points(manager, new_rings, first, prev_itr); + count = 0; + } +} + +template +void reassign_children_if_necessary(ring_ptr new_ring, + ring_ptr sibling_ring, + ring_manager& manager, + ring_vector& new_rings) { + auto& children = sibling_ring == nullptr ? manager.children : sibling_ring->children; + for (auto c : children) { + if (c == nullptr) { + continue; + } + if (std::find(new_rings.begin(), new_rings.end(), c) != new_rings.end()) { + continue; + } + if (poly2_contains_poly1(c, new_ring)) { + reassign_as_child(c, new_ring, manager); + } + } +} + +template +bool find_parent_in_tree(ring_ptr r, ring_ptr possible_parent, ring_manager& manager) { + // Before starting this we are assuming that possible_parent + // and r have opposite signs of their areas + + // First we must search all grandchildren + for (auto c : possible_parent->children) { + if (c == nullptr) { + continue; + } + for (auto gc : c->children) { + if (gc == nullptr) { + continue; + } + if (find_parent_in_tree(r, gc, manager)) { + return true; + } + } + } + + if (poly2_contains_poly1(r, possible_parent)) { + reassign_as_child(r, possible_parent, manager); + return true; + } + return false; +} + +template +void assign_new_ring_parents(ring_manager& manager, ring_ptr original_ring, ring_vector& new_rings) { + + // First lets remove any rings that have zero area + // or have no points + new_rings.erase(std::remove_if(new_rings.begin(), new_rings.end(), + [](ring_ptr const& r) { + if (r->points == nullptr) { + return true; + } + return value_is_zero(r->area()); + }), + new_rings.end()); + + if (new_rings.empty()) { + // No new rings created simply return; + return; + } + + // We should not have to re-assign the parent of the original ring + // because we always maintained the largest ring during splitting + // on repeated points. + + double original_ring_area = original_ring->area(); + bool original_positive = original_ring_area > 0.0; + + // If there is only one new ring the logic is very simple and we + // do not have to check which ring contains, we only need to compare + // the areas of the original ring and that of the new ring. + if (new_rings.size() == 1) { + double new_ring_area = new_rings.front()->area(); + bool new_positive = new_ring_area > 0.0; + if (original_positive == new_positive) { + // The rings should be siblings + assign_as_child(new_rings.front(), original_ring->parent, manager); + reassign_children_if_necessary(new_rings.front(), original_ring, manager, new_rings); + } else { + // The new ring is a child of original ring + // Check the + assign_as_child(new_rings.front(), original_ring, manager); + reassign_children_if_necessary(new_rings.front(), original_ring->parent, manager, new_rings); + } + return; + } + + // Now we want to sort rings from the largest in absolute area to the smallest + // as we will assign the rings with the largest areas first + std::stable_sort(new_rings.begin(), new_rings.end(), [](ring_ptr const& r1, ring_ptr const& r2) { + return std::fabs(r1->area()) > std::fabs(r2->area()); + }); + + for (auto r_itr = new_rings.begin(); r_itr != new_rings.end(); ++r_itr) { + double new_ring_area = (*r_itr)->area(); + bool new_positive = new_ring_area > 0.0; + bool same_orientation = new_positive == original_positive; + bool found = false; + // First lets check the trees of any new_rings that might have + // been assigned as siblings to the original ring. + for (auto s_itr = new_rings.begin(); s_itr != r_itr; ++s_itr) { + if ((*s_itr)->parent != original_ring->parent) { + continue; + } + if (same_orientation) { + for (auto s_child : (*s_itr)->children) { + if (s_child == nullptr) { + continue; + } + if (find_parent_in_tree(*r_itr, s_child, manager)) { + reassign_children_if_necessary(*r_itr, original_ring, manager, new_rings); + found = true; + break; + } + } + } else { + if (find_parent_in_tree(*r_itr, *s_itr, manager)) { + reassign_children_if_necessary(*r_itr, original_ring->parent, manager, new_rings); + found = true; + } + } + if (found) { + break; + } + } + + if (found) { + continue; + } + + // Next lets check the tree of the original_ring + if (same_orientation) { + for (auto o_child : original_ring->children) { + if (o_child == nullptr) { + continue; + } + if (find_parent_in_tree(*r_itr, o_child, manager)) { + reassign_children_if_necessary(*r_itr, original_ring, manager, new_rings); + found = true; + break; + } + } + if (!found) { + // If we didn't find any parent and the same orientation + // then it must be a sibling of the original ring + assign_as_child(*r_itr, original_ring->parent, manager); + reassign_children_if_necessary(*r_itr, original_ring, manager, new_rings); + } + } else { + if (find_parent_in_tree(*r_itr, original_ring, manager)) { + reassign_children_if_necessary(*r_itr, original_ring->parent, manager, new_rings); + } else { + throw std::runtime_error("Unable to find a proper parent ring"); + } + } + } +} + +template +bool correct_ring_self_intersections(ring_manager& manager, ring_ptr r, bool correct_tree) { + + if (r->corrected || !r->points) { + return false; + } + + ring_vector new_rings; + + find_and_correct_repeated_points(r, manager, new_rings); + + if (correct_tree) { + assign_new_ring_parents(manager, r, new_rings); + } + + r->corrected = true; + return true; +} + +template +void process_single_intersection(std::unordered_multimap, point_ptr_pair>& connection_map, + point_ptr op_j, + point_ptr op_k, + ring_manager& manager) { + ring_ptr ring_j = op_j->ring; + ring_ptr ring_k = op_k->ring; + if (ring_j == ring_k) { + return; + } + + if (!ring_j->is_hole() && !ring_k->is_hole()) { + // Both are not holes, return nothing to do. + return; + } + + ring_ptr ring_origin; + ring_ptr ring_search; + ring_ptr ring_parent; + point_ptr op_origin_1; + point_ptr op_origin_2; + if (!ring_j->is_hole()) { + ring_origin = ring_j; + ring_parent = ring_origin; + ring_search = ring_k; + op_origin_1 = op_j; + op_origin_2 = op_k; + } else if (!ring_k->is_hole()) { + ring_origin = ring_k; + ring_parent = ring_origin; + ring_search = ring_j; + op_origin_1 = op_k; + op_origin_2 = op_j; + + } else { + // both are holes + // Order doesn't matter + ring_origin = ring_j; + ring_parent = ring_origin->parent; + ring_search = ring_k; + op_origin_1 = op_j; + op_origin_2 = op_k; + } + if (ring_parent != ring_search->parent) { + // The two holes do not have the same parent, do not add them + // simply return! + return; + } + bool found = false; + std::list, point_ptr_pair>> iList; + { + auto range = connection_map.equal_range(ring_search); + // Check for direct connection + for (auto& it = range.first; it != range.second;) { + if (!it->second.op1->ring) { + it = connection_map.erase(it); + continue; + } + if (!it->second.op2->ring) { + it = connection_map.erase(it); + continue; + } + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring2 == ring_origin) { + found = true; + if (*op_origin_1 != *(it->second.op2)) { + iList.emplace_back(ring_search, it->second); + break; + } + } + ++it; + } + } + if (iList.empty()) { + auto range = connection_map.equal_range(ring_search); + std::set> visited; + visited.insert(ring_search); + // Check for connection through chain of other intersections + for (auto& it = range.first; it != range.second && it != connection_map.end() && it->first == ring_search; + ++it) { + ring_ptr it_ring = it->second.op2->ring; + if (it_ring != ring_search && *op_origin_2 != *it->second.op2 && it_ring != nullptr && + (ring_parent == it_ring || ring_parent == it_ring->parent) && !value_is_zero(it_ring->area()) && + find_intersect_loop(connection_map, iList, ring_parent, ring_origin, it_ring, visited, op_origin_2, + it->second.op2, manager)) { + found = true; + iList.emplace_front(ring_search, it->second); + break; + } + } + } + if (!found) { + point_ptr_pair intPt_origin = { op_origin_1, op_origin_2 }; + point_ptr_pair intPt_search = { op_origin_2, op_origin_1 }; + connection_map.emplace(ring_origin, std::move(intPt_origin)); + connection_map.emplace(ring_search, std::move(intPt_search)); + return; + } + + if (iList.empty()) { + // The situation where both origin and search are holes might have a missing + // search condition, we must check if a new pair must be added. + bool missing = true; + auto rng = connection_map.equal_range(ring_origin); + // Check for direct connection + for (auto& it = rng.first; it != rng.second; ++it) { + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring2 == ring_search) { + missing = false; + } + } + if (missing) { + point_ptr_pair intPt_origin = { op_origin_1, op_origin_2 }; + connection_map.emplace(ring_origin, std::move(intPt_origin)); + } + return; + } + if (ring_origin->is_hole()) { + for (auto& iRing : iList) { + ring_ptr ring_itr = iRing.first; + if (!ring_itr->is_hole()) { + // Make the hole the origin! + point_ptr op1 = op_origin_1; + op_origin_1 = iRing.second.op1; + iRing.second.op1 = op1; + point_ptr op2 = op_origin_2; + op_origin_2 = iRing.second.op2; + iRing.second.op2 = op2; + iRing.first = ring_origin; + ring_origin = ring_itr; + ring_parent = ring_origin; + break; + } + } + } + bool origin_is_hole = ring_origin->is_hole(); + + // Switch + point_ptr op_origin_1_next = op_origin_1->next; + point_ptr op_origin_2_next = op_origin_2->next; + op_origin_1->next = op_origin_2_next; + op_origin_2->next = op_origin_1_next; + op_origin_1_next->prev = op_origin_2; + op_origin_2_next->prev = op_origin_1; + + for (auto& iRing : iList) { + point_ptr op_search_1 = iRing.second.op1; + point_ptr op_search_2 = iRing.second.op2; + point_ptr op_search_1_next = op_search_1->next; + point_ptr op_search_2_next = op_search_2->next; + op_search_1->next = op_search_2_next; + op_search_2->next = op_search_1_next; + op_search_1_next->prev = op_search_2; + op_search_2_next->prev = op_search_1; + } + + ring_ptr ring_new = create_new_ring(manager); + ring_origin->corrected = false; + std::size_t size_1 = 0; + std::size_t size_2 = 0; + mapbox::geometry::box box1({ 0, 0 }, { 0, 0 }); + mapbox::geometry::box box2({ 0, 0 }, { 0, 0 }); + double area_1 = area_from_point(op_origin_1, size_1, box1); + double area_2 = area_from_point(op_origin_2, size_2, box2); + if (origin_is_hole && ((area_1 < 0.0))) { + ring_origin->points = op_origin_1; + ring_origin->set_stats(area_1, size_1, box1); + ring_new->points = op_origin_2; + ring_new->set_stats(area_2, size_2, box2); + } else { + ring_origin->points = op_origin_2; + ring_origin->set_stats(area_2, size_2, box2); + ring_new->points = op_origin_1; + ring_new->set_stats(area_1, size_1, box1); + } + + update_points_ring(ring_origin); + update_points_ring(ring_new); + + ring_origin->bottom_point = nullptr; + + for (auto& iRing : iList) { + ring_ptr ring_itr = iRing.first; + ring_itr->bottom_point = nullptr; + if (origin_is_hole) { + ring1_replaces_ring2(ring_origin, ring_itr, manager); + } else { + ring1_replaces_ring2(ring_origin->parent, ring_itr, manager); + } + } + if (origin_is_hole) { + assign_as_child(ring_new, ring_origin, manager); + // The parent ring in this situation might need to give up children + // to the new ring. + for (auto c : ring_parent->children) { + if (c == nullptr) { + continue; + } + if (poly2_contains_poly1(c, ring_new)) { + reassign_as_child(c, ring_new, manager); + } + } + } else { + // The new ring and the origin ring need to be siblings + // however some children ring from the ring origin might + // need to be re-assigned to the new ring + assign_as_sibling(ring_new, ring_origin, manager); + for (auto c : ring_origin->children) { + if (c == nullptr) { + continue; + } + if (poly2_contains_poly1(c, ring_new)) { + reassign_as_child(c, ring_new, manager); + } + } + } + + std::list, point_ptr_pair>> move_list; + + for (auto& iRing : iList) { + auto range_itr = connection_map.equal_range(iRing.first); + if (range_itr.first != range_itr.second) { + for (auto& it = range_itr.first; it != range_itr.second; ++it) { + ring_ptr it_ring = it->second.op1->ring; + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring == nullptr || it_ring2 == nullptr || it_ring == it_ring2) { + continue; + } + if (it_ring->is_hole() || it_ring2->is_hole()) { + move_list.emplace_back(it_ring, it->second); + } + } + connection_map.erase(iRing.first); + } + } + + auto range_itr = connection_map.equal_range(ring_origin); + for (auto& it = range_itr.first; it != range_itr.second;) { + ring_ptr it_ring = it->second.op1->ring; + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring == nullptr || it_ring2 == nullptr || it_ring == it_ring2) { + it = connection_map.erase(it); + continue; + } + if (it_ring != ring_origin) { + if (it_ring->is_hole() || it_ring2->is_hole()) { + move_list.emplace_back(it_ring, it->second); + } + it = connection_map.erase(it); + } else { + if (it_ring->is_hole() || it_ring2->is_hole()) { + ++it; + } else { + it = connection_map.erase(it); + } + } + } + + if (!move_list.empty()) { + connection_map.insert(move_list.begin(), move_list.end()); + } + + return; +} + +template +void correct_chained_repeats(ring_manager& manager, + std::unordered_multimap, point_ptr_pair>& connection_map, + point_vector_itr const& begin, + point_vector_itr const& end) { + for (auto itr1 = begin; itr1 != end; ++itr1) { + if ((*itr1)->ring == nullptr) { + continue; + } + for (auto itr2 = std::next(itr1); itr2 != end; ++itr2) { + if ((*itr2)->ring == nullptr) { + continue; + } + process_single_intersection(connection_map, *itr1, *itr2, manager); + } + } +} + +template +void correct_chained_rings(ring_manager& manager) { + + if (manager.all_points.size() < 2) { + return; + } + // Setup connection map which is a map of rings and their + // connection point pairs with other rings. + std::unordered_multimap, point_ptr_pair> connection_map; + connection_map.reserve(manager.rings.size()); + + // Now lets find and process any points + // that overlap -- we should have solved + // all situations where these points + // would be self intersections of a ring with + // earlier processing so this should just be + // points where different rings are touching. + std::size_t count = 0; + auto prev_itr = manager.all_points.begin(); + auto itr = std::next(prev_itr); + while (itr != manager.all_points.end()) { + if (*(*prev_itr) == *(*(itr))) { + ++count; + ++prev_itr; + ++itr; + if (itr != manager.all_points.end()) { + continue; + } else { + ++prev_itr; + } + } else { + ++prev_itr; + ++itr; + } + if (count == 0) { + continue; + } + auto first = prev_itr; + std::advance(first, -(static_cast(count) + 1)); + correct_chained_repeats(manager, connection_map, first, prev_itr); + count = 0; + } +} + +template +ring_vector sort_rings_largest_to_smallest(ring_manager& manager) { + ring_vector sorted_rings; + sorted_rings.reserve(manager.rings.size()); + for (auto& r : manager.rings) { + sorted_rings.push_back(&r); + } + std::stable_sort(sorted_rings.begin(), sorted_rings.end(), [](ring_ptr const& r1, ring_ptr const& r2) { + if (!r1->points || !r2->points) { + return r1->points != nullptr; + } + return std::fabs(r1->area()) > std::fabs(r2->area()); + }); + return sorted_rings; +} + +template +ring_vector sort_rings_smallest_to_largest(ring_manager& manager) { + ring_vector sorted_rings; + sorted_rings.reserve(manager.rings.size()); + for (auto& r : manager.rings) { + sorted_rings.push_back(&r); + } + std::stable_sort(sorted_rings.begin(), sorted_rings.end(), [](ring_ptr const& r1, ring_ptr const& r2) { + if (!r1->points || !r2->points) { + return r1->points != nullptr; + } + return std::fabs(r1->area()) < std::fabs(r2->area()); + }); + return sorted_rings; +} + +template +struct collinear_path { + // Collinear edges are in opposite directions + // such that start_1 is at the same position + // of end_2 and vise versa. Start to end is + // always forward (next) on path. + point_ptr start_1; + point_ptr end_1; + point_ptr start_2; + point_ptr end_2; +}; + +template +struct collinear_result { + point_ptr pt1; + point_ptr pt2; +}; + +template +collinear_result fix_collinear_path(collinear_path& path) { + // Left and right are just the opposite ends of the + // collinear path, they may not be actually left + // and right of each other. + // The left end is start_1 and end_2 + // The right end is start_2 and end_1 + + // NOTE spike detection is checking that the + // pointers are the same values, not position! + // additionally duplicate points we will treat + // if they are a spike left. + bool spike_left = (path.start_1 == path.end_2); + bool spike_right = (path.start_2 == path.end_1); + + if (spike_left && spike_right) { + // If both ends are spikes we should simply + // delete all the points. (they should be in a loop) + point_ptr itr = path.start_1; + while (itr != nullptr) { + itr->prev->next = nullptr; + itr->prev = nullptr; + itr->ring = nullptr; + itr = itr->next; + } + return { nullptr, nullptr }; + } else if (spike_left) { + point_ptr prev = path.start_2->prev; + point_ptr itr = path.start_2; + while (itr != path.end_1) { + itr->prev->next = nullptr; + itr->prev = nullptr; + itr->ring = nullptr; + itr = itr->next; + } + prev->next = path.end_1; + path.end_1->prev = prev; + return { path.end_1, nullptr }; + } else if (spike_right) { + point_ptr prev = path.start_1->prev; + point_ptr itr = path.start_1; + while (itr != path.end_2) { + itr->prev->next = nullptr; + itr->prev = nullptr; + itr->ring = nullptr; + itr = itr->next; + } + prev->next = path.end_2; + path.end_2->prev = prev; + return { path.end_2, nullptr }; + } else { + point_ptr prev_1 = path.start_1->prev; + point_ptr prev_2 = path.start_2->prev; + point_ptr itr = path.start_1; + do { + itr->prev->next = nullptr; + itr->prev = nullptr; + itr->ring = nullptr; + itr = itr->next; + } while (itr != path.end_1 && itr != nullptr); + itr = path.start_2; + do { + itr->prev->next = nullptr; + itr->prev = nullptr; + itr->ring = nullptr; + itr = itr->next; + } while (itr != path.end_2 && itr != nullptr); + if (path.start_1 == path.end_1 && path.start_2 == path.end_2) { + return { nullptr, nullptr }; + } else if (path.start_1 == path.end_1) { + prev_2->next = path.end_2; + path.end_2->prev = prev_2; + return { path.end_2, nullptr }; + } else if (path.start_2 == path.end_2) { + prev_1->next = path.end_1; + path.end_1->prev = prev_1; + return { path.end_1, nullptr }; + } else { + prev_1->next = path.end_2; + path.end_2->prev = prev_1; + prev_2->next = path.end_1; + path.end_1->prev = prev_2; + return { path.end_1, path.end_2 }; + } + } +} + +template +collinear_path find_start_and_end_of_collinear_edges(point_ptr pt_a, point_ptr pt_b) { + // Search backward on A, forwards on B first + bool same_ring = pt_a->ring == pt_b->ring; + point_ptr back = pt_a; + point_ptr forward = pt_b; + bool first = true; + do { + while (*(back->prev) == *back && back != forward) { + back = back->prev; + if (back == pt_a) { + break; + } + } + if (back == forward) { + back = back->prev; + forward = forward->next; + break; + } + while (*(forward->next) == *forward && back != forward) { + forward = forward->next; + if (forward == pt_b) { + break; + } + } + if (!first && (back == pt_a || forward == pt_b)) { + break; + } + if (back == forward) { + back = back->prev; + forward = forward->next; + break; + } + back = back->prev; + forward = forward->next; + first = false; + } while (*back == *forward); + point_ptr start_a = back->next; + // If there are repeated points at the diverge we want to select + // only the first of those repeated points. + while (!same_ring && *start_a == *start_a->next && start_a != pt_a) { + start_a = start_a->next; + } + point_ptr end_b = forward->prev; + while (!same_ring && *end_b == *end_b->prev && end_b != pt_b) { + end_b = end_b->prev; + } + // Search backward on B, forwards on A next + back = pt_b; + forward = pt_a; + first = true; + do { + while (*(back->prev) == *back && back != forward) { + back = back->prev; + if (back == pt_b) { + break; + } + } + if (back == forward) { + back = back->prev; + forward = forward->next; + break; + } + while (*(forward->next) == *forward && back != forward) { + forward = forward->next; + if (forward == pt_a) { + break; + } + } + if (!first && (back == pt_b || forward == pt_a)) { + break; + } + if (back == forward || (!first && (back == end_b || forward == start_a))) { + back = back->prev; + forward = forward->next; + break; + } + back = back->prev; + forward = forward->next; + first = false; + } while (*back == *forward); + point_ptr start_b = back->next; + while (!same_ring && *start_b == *start_b->next && start_b != pt_b) { + start_b = start_b->next; + } + point_ptr end_a = forward->prev; + while (!same_ring && *end_a == *end_a->prev && end_a != pt_a) { + end_a = end_a->prev; + } + return { start_a, end_a, start_b, end_b }; +} + +template +bool has_collinear_edge(point_ptr pt_a, point_ptr pt_b) { + // It is assumed pt_a and pt_b are at the same location. + return (*pt_a->next == *pt_b->prev || *pt_b->next == *pt_a->prev); +} + +template +void process_collinear_edges_same_ring(point_ptr pt_a, point_ptr pt_b, ring_manager& manager) { + ring_ptr original_ring = pt_a->ring; + // As they are the same ring that are forming a collinear edge + // we should expect the creation of two different rings. + auto path = find_start_and_end_of_collinear_edges(pt_a, pt_b); + auto results = fix_collinear_path(path); + if (results.pt1 == nullptr) { + // If pt1 is a nullptr, both values + // are nullptrs. This mean the ring was + // removed during this processing. + remove_ring(original_ring, manager, false); + } else if (results.pt2 == nullptr) { + // If a single point is only returned, we simply removed a spike. + // In this case, we don't need to worry about parent or children + // and we simply need to set the points and clear the area + original_ring->points = results.pt1; + original_ring->recalculate_stats(); + } else { + // If we have two seperate points, the ring has split into + // two different rings. + ring_ptr ring_new = create_new_ring(manager); + ring_new->points = results.pt2; + ring_new->recalculate_stats(); + update_points_ring(ring_new); + original_ring->points = results.pt1; + original_ring->recalculate_stats(); + } +} + +template +void process_collinear_edges_different_rings(point_ptr pt_a, point_ptr pt_b, ring_manager& manager) { + ring_ptr ring_a = pt_a->ring; + ring_ptr ring_b = pt_b->ring; + bool ring_a_larger = std::fabs(ring_a->area()) > std::fabs(ring_b->area()); + auto path = find_start_and_end_of_collinear_edges(pt_a, pt_b); + // This should result in two rings becoming one. + auto results = fix_collinear_path(path); + if (results.pt1 == nullptr) { + remove_ring(ring_a, manager, false); + remove_ring(ring_b, manager, false); + return; + } + // Rings should merge into a single ring of the same orientation. + // Therefore, we we will need to replace one ring with the other + ring_ptr merged_ring = ring_a_larger ? ring_a : ring_b; + ring_ptr deleted_ring = ring_a_larger ? ring_b : ring_a; + + merged_ring->points = results.pt1; + update_points_ring(merged_ring); + merged_ring->recalculate_stats(); + if (merged_ring->size() < 3) { + remove_ring_and_points(merged_ring, manager, false); + } + remove_ring(deleted_ring, manager, false); +} + +template +bool remove_duplicate_points(point_ptr pt_a, point_ptr pt_b, ring_manager& manager) { + if (pt_a->ring == pt_b->ring) { + if (pt_a->next == pt_b) { + pt_a->next = pt_b->next; + pt_a->next->prev = pt_a; + pt_b->next = nullptr; + pt_b->prev = nullptr; + pt_b->ring = nullptr; + if (pt_a->ring->points == pt_b) { + pt_a->ring->points = pt_a; + } + return true; + } else if (pt_b->next == pt_a) { + pt_a->prev = pt_b->prev; + pt_a->prev->next = pt_a; + pt_b->next = nullptr; + pt_b->prev = nullptr; + pt_b->ring = nullptr; + if (pt_a->ring->points == pt_b) { + pt_a->ring->points = pt_a; + } + return true; + } + } + while (*pt_a->next == *pt_a && pt_a->next != pt_a) { + point_ptr remove = pt_a->next; + pt_a->next = remove->next; + pt_a->next->prev = pt_a; + remove->next = nullptr; + remove->prev = nullptr; + remove->ring = nullptr; + if (pt_a->ring->points == remove) { + pt_a->ring->points = pt_a; + } + } + while (*pt_a->prev == *pt_a && pt_a->prev != pt_a) { + point_ptr remove = pt_a->prev; + pt_a->prev = remove->prev; + pt_a->prev->next = pt_a; + remove->next = nullptr; + remove->prev = nullptr; + remove->ring = nullptr; + if (pt_a->ring->points == remove) { + pt_a->ring->points = pt_a; + } + } + if (pt_a->next == pt_a) { + remove_ring_and_points(pt_a->ring, manager, false); + return true; + } + if (pt_b->ring == nullptr) { + return true; + } + while (*pt_b->next == *pt_b && pt_b->next != pt_b) { + point_ptr remove = pt_b->next; + pt_b->next = remove->next; + pt_b->next->prev = pt_b; + remove->next = nullptr; + remove->prev = nullptr; + remove->ring = nullptr; + if (pt_b->ring->points == remove) { + pt_b->ring->points = pt_b; + } + } + while (*pt_b->prev == *pt_b && pt_b->prev != pt_b) { + point_ptr remove = pt_b->prev; + pt_b->prev = remove->prev; + pt_b->prev->next = pt_b; + remove->next = nullptr; + remove->prev = nullptr; + remove->ring = nullptr; + if (pt_b->ring->points == remove) { + pt_b->ring->points = pt_b; + } + } + if (pt_b->next == pt_b) { + remove_ring_and_points(pt_b->ring, manager, false); + return true; + } + if (pt_a->ring == nullptr) { + return true; + } + return false; +} + +template +bool process_collinear_edges(point_ptr pt_a, point_ptr pt_b, ring_manager& manager) { + // Neither point assigned to a ring (deleted points) + if (!pt_a->ring || !pt_b->ring) { + return false; + } + + if (remove_duplicate_points(pt_a, pt_b, manager)) { + return true; + } + + if (!has_collinear_edge(pt_a, pt_b)) { + if (pt_a->ring == pt_b->ring) { + correct_self_intersection(pt_a, pt_b, manager); + return true; + } + return false; + } + + if (pt_a->ring == pt_b->ring) { + process_collinear_edges_same_ring(pt_a, pt_b, manager); + } else { + process_collinear_edges_different_rings(pt_a, pt_b, manager); + } + return true; +} + +template +void correct_collinear_repeats(ring_manager& manager, + point_vector_itr const& begin, + point_vector_itr const& end) { + for (auto itr1 = begin; itr1 != end; ++itr1) { + if ((*itr1)->ring == nullptr) { + continue; + } + for (auto itr2 = begin; itr2 != end;) { + if ((*itr1)->ring == nullptr) { + break; + } + if ((*itr2)->ring == nullptr || *itr2 == *itr1) { + ++itr2; + continue; + } + if (process_collinear_edges(*itr1, *itr2, manager)) { + itr2 = begin; + } else { + ++itr2; + } + } + } +} + +template +void correct_collinear_edges(ring_manager& manager) { + + if (manager.all_points.size() < 2) { + return; + } + std::size_t count = 0; + auto prev_itr = manager.all_points.begin(); + auto itr = std::next(prev_itr); + while (itr != manager.all_points.end()) { + if (*(*prev_itr) == *(*(itr))) { + ++count; + ++prev_itr; + ++itr; + if (itr != manager.all_points.end()) { + continue; + } else { + ++prev_itr; + } + } else { + ++prev_itr; + ++itr; + } + if (count == 0) { + continue; + } + auto first = prev_itr; + std::advance(first, -(static_cast(count) + 1)); + correct_collinear_repeats(manager, first, prev_itr); + count = 0; + } +} + +template +void correct_tree(ring_manager& manager) { + + // It is possible that vatti resulted in some parent child + // relationships that are not quite right, therefore, we + // need to just rebuild the entire tree of rings + + // First lets process the rings in order of size from largest + // area to smallest, we know right away that no smaller ring could ever + // contain a larger ring so we can use this to our advantage + // as we iterate over the rings. + using rev_itr = typename ring_vector::reverse_iterator; + ring_vector sorted_rings = sort_rings_largest_to_smallest(manager); + for (auto itr = sorted_rings.begin(); itr != sorted_rings.end(); ++itr) { + if ((*itr)->points == nullptr) { + continue; + } + if ((*itr)->size() < 3 || value_is_zero((*itr)->area())) { + remove_ring_and_points(*itr, manager, false); + continue; + } + (*itr)->corrected = true; + bool found = false; + // Search in reverse from the current iterator back to the begining + // to see if any of those rings might be its parent. + for (auto r_itr = rev_itr(itr); r_itr != sorted_rings.rend(); ++r_itr) { + // If orientations are not different, this can't be its parent. + if ((*r_itr)->is_hole() == (*itr)->is_hole()) { + continue; + } + if (poly2_contains_poly1(*itr, *r_itr)) { + reassign_as_child(*itr, *r_itr, manager); + found = true; + break; + } + } + if (!found) { + if ((*itr)->is_hole()) { + throw std::runtime_error("Could not properly place hole to a parent."); + } else { + // Assign to base of tree by passing nullptr + reassign_as_child(*itr, static_cast>(nullptr), manager); + } + } + } +} + +template +bool correct_self_intersections(ring_manager& manager, bool correct_tree) { + bool fixed_intersections = false; + auto sorted_rings = sort_rings_smallest_to_largest(manager); + for (auto const& r : sorted_rings) { + if (correct_ring_self_intersections(manager, r, correct_tree)) { + fixed_intersections = true; + } + } + return fixed_intersections; +} + +template +void correct_topology(ring_manager& manager) { + + // Sort all the points, this will be used for the locating of chained rings + // and the collinear edges and only needs to be done once. + std::stable_sort(manager.all_points.begin(), manager.all_points.end(), point_ptr_cmp()); + + // Initially the orientations of the rings + // could be incorrect, we need to adjust them + correct_orientations(manager); + + // We should only have to fix collinear edges once. + // During this we also correct self intersections + correct_collinear_edges(manager); + + correct_self_intersections(manager, false); + + correct_tree(manager); + + bool fixed_intersections = true; + while (fixed_intersections) { + correct_chained_rings(manager); + fixed_intersections = correct_self_intersections(manager, true); + } +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/util.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/util.hpp new file mode 100644 index 000000000..f8acdb476 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/util.hpp @@ -0,0 +1,89 @@ +#pragma once + +#include + +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +double area(mapbox::geometry::linear_ring const& poly) { + std::size_t size = poly.size(); + if (size < 3) { + return 0.0; + } + + double a = 0.0; + auto itr = poly.begin(); + auto itr_prev = poly.end(); + --itr_prev; + a += static_cast(itr_prev->x + itr->x) * static_cast(itr_prev->y - itr->y); + ++itr; + itr_prev = poly.begin(); + for (; itr != poly.end(); ++itr, ++itr_prev) { + a += static_cast(itr_prev->x + itr->x) * static_cast(itr_prev->y - itr->y); + } + return -a * 0.5; +} + +inline bool value_is_zero(double val) { + return std::fabs(val) < (5.0 * std::numeric_limits::epsilon()); +} + +inline bool values_are_equal(double x, double y) { + return value_is_zero(x - y); +} + +inline bool greater_than_or_equal(double x, double y) { + return x > y || values_are_equal(x, y); +} + +template +bool slopes_equal(mapbox::geometry::point const& pt1, + mapbox::geometry::point const& pt2, + mapbox::geometry::point const& pt3) { + return static_cast(pt1.y - pt2.y) * static_cast(pt2.x - pt3.x) == + static_cast(pt1.x - pt2.x) * static_cast(pt2.y - pt3.y); +} + +template +bool slopes_equal(mapbox::geometry::wagyu::point const& pt1, + mapbox::geometry::wagyu::point const& pt2, + mapbox::geometry::point const& pt3) { + return static_cast(pt1.y - pt2.y) * static_cast(pt2.x - pt3.x) == + static_cast(pt1.x - pt2.x) * static_cast(pt2.y - pt3.y); +} + +template +bool slopes_equal(mapbox::geometry::wagyu::point const& pt1, + mapbox::geometry::wagyu::point const& pt2, + mapbox::geometry::wagyu::point const& pt3) { + return static_cast(pt1.y - pt2.y) * static_cast(pt2.x - pt3.x) == + static_cast(pt1.x - pt2.x) * static_cast(pt2.y - pt3.y); +} + +template +bool slopes_equal(mapbox::geometry::point const& pt1, + mapbox::geometry::point const& pt2, + mapbox::geometry::point const& pt3, + mapbox::geometry::point const& pt4) { + return static_cast(pt1.y - pt2.y) * static_cast(pt3.x - pt4.x) == + static_cast(pt1.x - pt2.x) * static_cast(pt3.y - pt4.y); +} + +template +inline T wround(double value) { + return static_cast(::llround(value)); +} + +template <> +inline std::int64_t wround(double value) { + return ::llround(value); +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/vatti.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/vatti.hpp new file mode 100644 index 000000000..6f675ea52 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/vatti.hpp @@ -0,0 +1,64 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +void execute_vatti(local_minimum_list& minima_list, + ring_manager& manager, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + active_bound_list active_bounds; + scanbeam_list scanbeam; + T scanline_y = std::numeric_limits::max(); + + local_minimum_ptr_list minima_sorted; + minima_sorted.reserve(minima_list.size()); + for (auto& lm : minima_list) { + minima_sorted.push_back(&lm); + } + std::stable_sort(minima_sorted.begin(), minima_sorted.end(), local_minimum_sorter()); + local_minimum_ptr_list_itr current_lm = minima_sorted.begin(); + // std::clog << output_all_edges(minima_sorted) << std::endl; + + setup_scanbeam(minima_list, scanbeam); + manager.current_hp_itr = manager.hot_pixels.begin(); + + while (pop_from_scanbeam(scanline_y, scanbeam) || current_lm != minima_sorted.end()) { + + process_intersections(scanline_y, active_bounds, cliptype, subject_fill_type, clip_fill_type, manager); + + update_current_hp_itr(scanline_y, manager); + + // First we process bounds that has already been added to the active bound list -- + // if the active bound list is empty local minima that are at this scanline_y and + // have a horizontal edge at the local minima will be processed + process_edges_at_top_of_scanbeam(scanline_y, active_bounds, scanbeam, minima_sorted, current_lm, manager, + cliptype, subject_fill_type, clip_fill_type); + + // Next we will add local minima bounds to the active bounds list that are on the local + // minima queue at + // this current scanline_y + insert_local_minima_into_ABL(scanline_y, minima_sorted, current_lm, active_bounds, manager, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + } +} +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry/wagyu/wagyu.hpp b/deps/wagyu/include/mapbox/geometry/wagyu/wagyu.hpp new file mode 100644 index 000000000..440141810 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry/wagyu/wagyu.hpp @@ -0,0 +1,145 @@ +#pragma once + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#define WAGYU_MAJOR_VERSION 0 +#define WAGYU_MINOR_VERSION 4 +#define WAGYU_PATCH_VERSION 3 + +#define WAGYU_VERSION (WAGYU_MAJOR_VERSION * 100000) + (WAGYU_MINOR_VERSION * 100) + (WAGYU_PATCH_VERSION) + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +class wagyu { +private: + local_minimum_list minima_list; + bool reverse_output; + + wagyu(wagyu const&) = delete; + wagyu& operator=(wagyu const&) = delete; + +public: + wagyu() : minima_list(), reverse_output(false) { + } + + ~wagyu() { + clear(); + } + + template + bool add_ring(mapbox::geometry::linear_ring const& pg, polygon_type p_type = polygon_type_subject) { + return add_linear_ring(pg, minima_list, p_type); + } + + template + bool add_polygon(mapbox::geometry::polygon const& ppg, polygon_type p_type = polygon_type_subject) { + bool result = false; + for (auto const& r : ppg) { + if (add_ring(r, p_type)) { + result = true; + } + } + return result; + } + + void reverse_rings(bool value) { + reverse_output = value; + } + + void clear() { + minima_list.clear(); + } + + mapbox::geometry::box get_bounds() { + mapbox::geometry::point min = { 0, 0 }; + mapbox::geometry::point max = { 0, 0 }; + if (minima_list.empty()) { + return mapbox::geometry::box(min, max); + } + bool first_set = false; + for (auto const& lm : minima_list) { + if (!lm.left_bound.edges.empty()) { + if (!first_set) { + min = lm.left_bound.edges.front().top; + max = lm.left_bound.edges.back().bot; + first_set = true; + } else { + min.y = std::min(min.y, lm.left_bound.edges.front().top.y); + max.y = std::max(max.y, lm.left_bound.edges.back().bot.y); + max.x = std::max(max.x, lm.left_bound.edges.back().top.x); + min.x = std::min(min.x, lm.left_bound.edges.back().top.x); + } + for (auto const& e : lm.left_bound.edges) { + max.x = std::max(max.x, e.bot.x); + min.x = std::min(min.x, e.bot.x); + } + } + if (!lm.right_bound.edges.empty()) { + if (!first_set) { + min = lm.right_bound.edges.front().top; + max = lm.right_bound.edges.back().bot; + first_set = true; + } else { + min.y = std::min(min.y, lm.right_bound.edges.front().top.y); + max.y = std::max(max.y, lm.right_bound.edges.back().bot.y); + max.x = std::max(max.x, lm.right_bound.edges.back().top.x); + min.x = std::min(min.x, lm.right_bound.edges.back().top.x); + } + for (auto const& e : lm.right_bound.edges) { + max.x = std::max(max.x, e.bot.x); + min.x = std::min(min.x, e.bot.x); + } + } + } + return mapbox::geometry::box(min, max); + } + + template + bool execute(clip_type cliptype, + mapbox::geometry::multi_polygon& solution, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + if (minima_list.empty()) { + return false; + } + + ring_manager manager; + + interrupt_check(); // Check for interruptions + + build_hot_pixels(minima_list, manager); + + interrupt_check(); // Check for interruptions + + execute_vatti(minima_list, manager, cliptype, subject_fill_type, clip_fill_type); + + interrupt_check(); // Check for interruptions + + correct_topology(manager); + + build_result(solution, manager, reverse_output); + + return true; + } +}; +} // namespace wagyu +} // namespace geometry +} // namespace mapbox diff --git a/deps/wagyu/include/mapbox/geometry_io.hpp b/deps/wagyu/include/mapbox/geometry_io.hpp new file mode 100644 index 000000000..4d91ccc96 --- /dev/null +++ b/deps/wagyu/include/mapbox/geometry_io.hpp @@ -0,0 +1,98 @@ +#pragma once + +#include +#include + +#include +#include + +namespace mapbox { +namespace geometry { + +std::ostream& operator<<(std::ostream& os, const empty&) +{ + return os << "[]"; +} + +template +std::ostream& operator<<(std::ostream& os, const point& point) +{ + return os << "[" << point.x << "," << point.y << "]"; +} + +template class C, class... Args> +std::ostream& operator<<(std::ostream& os, const C& cont) +{ + os << "["; + for (auto it = cont.cbegin();;) + { + os << *it; + if (++it == cont.cend()) + { + break; + } + os << ","; + } + return os << "]"; +} + +template +std::ostream& operator<<(std::ostream& os, const line_string& geom) +{ + return os << static_cast::container_type>(geom); +} + +template +std::ostream& operator<<(std::ostream& os, const linear_ring& geom) +{ + return os << static_cast::container_type>(geom); +} + +template +std::ostream& operator<<(std::ostream& os, const polygon& geom) +{ + return os << static_cast::container_type>(geom); +} + +template +std::ostream& operator<<(std::ostream& os, const multi_point& geom) +{ + return os << static_cast::container_type>(geom); +} + +template +std::ostream& operator<<(std::ostream& os, const multi_line_string& geom) +{ + return os << static_cast::container_type>(geom); +} + +template +std::ostream& operator<<(std::ostream& os, const multi_polygon& geom) +{ + return os << static_cast::container_type>(geom); +} + +template +std::ostream& operator<<(std::ostream& os, const geometry& geom) +{ + geometry::visit(geom, [&](const auto& g) { os << g; }); + return os; +} + +template +std::ostream& operator<<(std::ostream& os, const geometry_collection& geom) +{ + return os << static_cast::container_type>(geom); +} + +} // namespace geometry + +namespace feature { + +std::ostream& operator<<(std::ostream& os, const null_value_t&) +{ + return os << "[]"; +} + +} // namespace feature +} // namespace mapbox diff --git a/deps/wagyu/lwgeom_wagyu.cpp b/deps/wagyu/lwgeom_wagyu.cpp new file mode 100644 index 000000000..c1c89e546 --- /dev/null +++ b/deps/wagyu/lwgeom_wagyu.cpp @@ -0,0 +1,296 @@ +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * + * PostGIS is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * PostGIS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with PostGIS. If not, see . + * + ********************************************************************** + * + * Copyright 2019 Raúl Marín + * + **********************************************************************/ + +#include "lwgeom_wagyu.h" + +#define USE_WAGYU_INTERRUPT +#include "mapbox/geometry/wagyu/wagyu.hpp" +#include "mapbox/geometry/wagyu/quick_clip.hpp" + +#include +#include +#include + +using namespace mapbox::geometry; + +using wagyu_coord_type = std::int32_t; +using wagyu_polygon = mapbox::geometry::polygon; +using wagyu_multipolygon = mapbox::geometry::multi_polygon; +using wagyu_linearring = mapbox::geometry::linear_ring; +using vwpolygon = std::vector; +using wagyu_point = mapbox::geometry::point; +using wagyu_box = mapbox::geometry::box; + +static wagyu_linearring +ptarray_to_wglinearring(const POINTARRAY *pa, const wagyu_box &tile) +{ + wagyu_linearring lr; + lr.reserve(pa->npoints); + + /* We calculate the bounding box of the point array */ + wagyu_coord_type min_x = std::numeric_limits::max(); + wagyu_coord_type max_x = std::numeric_limits::min(); + wagyu_coord_type min_y = min_x; + wagyu_coord_type max_y = max_x; + + for (std::uint32_t i = 0; i < pa->npoints; i++) + { + const POINT2D *p2d = getPoint2d_cp(pa, i); + wagyu_point wp(static_cast(p2d->x), static_cast(p2d->y)); + if (wp.x < min_x) + { + min_x = wp.x; + } + else if (wp.x > max_x) + { + max_x = wp.x; + } + + if (wp.y < min_y) + { + min_y = wp.y; + } + else if (wp.y > max_y) + { + max_y = wp.y; + } + + lr.push_back(std::move(wp)); + } + + /* Check how many sides of the calculated box are inside the tile */ + uint32_t sides_in = 0; + sides_in += min_x >= tile.min.x && min_x <= tile.max.x; + sides_in += max_x >= tile.min.x && max_x <= tile.max.x; + sides_in += min_y >= tile.min.y && min_y <= tile.max.y; + sides_in += max_y >= tile.min.y && max_y <= tile.max.y; + + /* With 0 sides in, the box it's either outside or covers the tile completely */ + if ((sides_in == 0) && (min_x > tile.max.x || max_x < tile.min.x || min_y > tile.max.y || max_y < tile.min.y)) + { + /* No overlapping: Return an empty linearring */ + return wagyu_linearring(); + } + + if (sides_in != 4) + { + /* Some edges need to be clipped */ + return mapbox::geometry::wagyu::quick_clip::quick_lr_clip(lr, tile); + } + + /* All points are inside the box */ + return lr; +} + +static vwpolygon +lwpoly_to_vwgpoly(const LWPOLY *geom, const wagyu_box &box) +{ + vwpolygon vp; + for (std::uint32_t i = 0; i < geom->nrings; i++) + { + wagyu_polygon pol; + pol.push_back(ptarray_to_wglinearring(geom->rings[i], box)); + + /* If it has an inner ring, push it too */ + i++; + if (i != geom->nrings) + { + pol.push_back(ptarray_to_wglinearring(geom->rings[i], box)); + } + + vp.push_back(pol); + } + + return vp; +} + +static vwpolygon +lwmpoly_to_vwgpoly(const LWMPOLY *geom, const wagyu_box &box) +{ + vwpolygon vp; + for (std::uint32_t i = 0; i < geom->ngeoms; i++) + { + vwpolygon vp_intermediate = lwpoly_to_vwgpoly(geom->geoms[i], box); + vp.insert(vp.end(), + std::make_move_iterator(vp_intermediate.begin()), + std::make_move_iterator(vp_intermediate.end())); + } + + return vp; +} + +/** + * Transforms a liblwgeom geometry (types POLYGONTYPE or MULTIPOLYGONTYPE) + * into an array of mapbox::geometry (polygon) + * A single lwgeom polygon can lead to N mapbox polygon as odd numbered rings + * are treated as new polygons (and even numbered treated as holes in the + * previous rings) + */ +static vwpolygon +lwgeom_to_vwgpoly(const LWGEOM *geom, const wagyu_box &box) +{ + switch (geom->type) + { + case POLYGONTYPE: + return lwpoly_to_vwgpoly(reinterpret_cast(geom), box); + case MULTIPOLYGONTYPE: + return lwmpoly_to_vwgpoly(reinterpret_cast(geom), box); + default: + return vwpolygon(); + } +} + +static POINTARRAY * +wglinearring_to_ptarray(const wagyu_polygon::linear_ring_type &lr) +{ + const uint32_t npoints = lr.size(); + POINTARRAY *pa = ptarray_construct(0, 0, npoints); + + for (uint32_t i = 0; i < npoints; i++) + { + const wagyu_polygon::linear_ring_type::point_type &p = lr[i]; + POINT4D point = {static_cast(p.x), static_cast(p.y), 0.0, 0.0}; + ptarray_set_point4d(pa, i, &point); + } + + return pa; +} + +static LWGEOM * +wgpoly_to_lwgeom(const wagyu_multipolygon::polygon_type &p) +{ + const uint32_t nrings = p.size(); + POINTARRAY **ppa = reinterpret_cast(lwalloc(sizeof(POINTARRAY *) * nrings)); + + for (uint32_t i = 0; i < nrings; i++) + { + ppa[i] = wglinearring_to_ptarray(p[i]); + } + + return reinterpret_cast(lwpoly_construct(0, NULL, nrings, ppa)); +} + +static LWGEOM * +wgmpoly_to_lwgeom(const wagyu_multipolygon &mp) +{ + const uint32_t ngeoms = mp.size(); + if (ngeoms == 1) + { + return wgpoly_to_lwgeom(mp[0]); + } + else + { + LWGEOM **geoms = reinterpret_cast(lwalloc(sizeof(LWGEOM *) * ngeoms)); + + for (uint32_t i = 0; i < ngeoms; i++) + { + geoms[i] = wgpoly_to_lwgeom(mp[i]); + } + + return reinterpret_cast(lwcollection_construct(MULTIPOLYGONTYPE, 0, NULL, ngeoms, geoms)); + } +} + +static LWGEOM * +__lwgeom_wagyu_clip_by_box(const LWGEOM *geom, const GBOX *bbox) +{ + if (!geom || !bbox) + { + return NULL; + } + + if (geom->type != POLYGONTYPE && geom->type != MULTIPOLYGONTYPE) + { + return NULL; + } + + if (lwgeom_is_empty(geom)) + { + LWGEOM *out = lwgeom_construct_empty(MULTIPOLYGONTYPE, geom->srid, 0, 0); + out->flags = geom->flags; + return out; + } + + wagyu_point box_min(std::min(bbox->xmin, bbox->xmax), std::min(bbox->ymin, bbox->ymax)); + wagyu_point box_max(std::max(bbox->xmin, bbox->xmax), std::max(bbox->ymin, bbox->ymax)); + const wagyu_box box(box_min, box_max); + + wagyu_multipolygon solution; + vwpolygon vpsubject = lwgeom_to_vwgpoly(geom, box); + if (vpsubject.size() == 0) + { + LWGEOM *out = lwgeom_construct_empty(MULTIPOLYGONTYPE, geom->srid, 0, 0); + out->flags = geom->flags; + return out; + } + + wagyu::wagyu clipper; + for (auto &poly : vpsubject) + { + for (auto &lr : poly) + { + if (!lr.empty()) + { + clipper.add_ring(lr, wagyu::polygon_type_subject); + } + } + } + + clipper.execute(wagyu::clip_type_union, solution, wagyu::fill_type_even_odd, wagyu::fill_type_even_odd); + + LWGEOM *g = wgmpoly_to_lwgeom(solution); + g->srid = geom->srid; + + return g; +} + +LWGEOM * +lwgeom_wagyu_clip_by_box(const LWGEOM *geom, const GBOX *bbox) +{ + mapbox::geometry::wagyu::interrupt_reset(); + try + { + return __lwgeom_wagyu_clip_by_box(geom, bbox); + } + catch (...) + { + return NULL; + } +} + +const char * +libwagyu_version() +{ + static char str[50] = {0}; + snprintf( + str, sizeof(str), "%d.%d.%d (Internal)", WAGYU_MAJOR_VERSION, WAGYU_MINOR_VERSION, WAGYU_PATCH_VERSION); + + return str; +} + +void +lwgeom_wagyu_interruptRequest() +{ + mapbox::geometry::wagyu::interrupt_request(); +} diff --git a/deps/wagyu/lwgeom_wagyu.h b/deps/wagyu/lwgeom_wagyu.h new file mode 100644 index 000000000..738ecee46 --- /dev/null +++ b/deps/wagyu/lwgeom_wagyu.h @@ -0,0 +1,66 @@ +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * + * PostGIS is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * PostGIS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with PostGIS. If not, see . + * + ********************************************************************** + * + * Copyright 2019 Raúl Marín + * + **********************************************************************/ + +#ifndef LWGEOM_WAGYU_H +#define LWGEOM_WAGYU_H + +#ifdef __cplusplus +extern "C" { +#endif + +#include + +/** + * Returns a geometry that represents the point set intersection of 2 polygons. + * The result is guaranteed to be valid according to OGC standard. + * This function only supports 2D and will DROP any extra dimensions. + * + * The return value will be either NULL (on error) or a new geometry with the + * same SRID as the first input. + * The ownership of it is given to the caller. + * + * @param geom - A geometry of either POLYGONTYPE or MULTIPOLYGONTYPE type. + * @param bbox - The GBOX that will be used to clip. + * @return - NULL on invalid input (NULL or wrong type), or when interrupted. + * - An empty MULTIPOLYGONTYPE the geometry is empty. + * - A pointer to a LWMPOLY otherwise. + */ +LWGEOM *lwgeom_wagyu_clip_by_box(const LWGEOM *geom, const GBOX *bbox); + +/** + * Returns a pointer to the static string representing the Wagyu release + * being used (Should not be freed) + */ +const char *libwagyu_version(); + +/** + * Requests wagyu to stop processing. In this case it will return NULL + */ +void lwgeom_wagyu_interruptRequest(); + +#ifdef __cplusplus +} +#endif + +#endif /* LWGEOM_WAGYU_H */ diff --git a/doc/installation.xml b/doc/installation.xml index 1f691c418..53d280a11 100644 --- a/doc/installation.xml +++ b/doc/installation.xml @@ -195,6 +195,7 @@ psql -d yourdatabase -f sfcgal_comments.sql https://trac.osgeo.org/postgis/ticket/4125. + @@ -240,6 +241,7 @@ psql -d yourdatabase -f sfcgal_comments.sql To enable ST_AsMVT protobuf-c library (for usage) and the protoc-c compiler (for building) are required. Also, pkg-config is required to verify the correct minimum version of protobuf-c. See protobuf-c. + To use Wagyu to validate MVT polygons faster, a c++11 compiler is required. It will use the CXX and CXXFLAGS and needs --with-wagyu to be passed during configure. @@ -596,6 +598,14 @@ tar -xvzf postgis-&last_release_version;.tar.gz + + --with-wagyu + + + When building with MVT support, by default Postgis will use GEOS to clip and validate MVT polygons. You can use Wagyu instead which is faster and guaranteed to produce correct values for this specific case. + + + diff --git a/doc/reference_management.xml b/doc/reference_management.xml index e7f112382..e28c74ea8 100644 --- a/doc/reference_management.xml +++ b/doc/reference_management.xml @@ -451,9 +451,10 @@ LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" TOPOLOGY SELECT PostGIS_Full_Version(); postgis_full_version ---------------------------------------------------------------------------------- -POSTGIS="2.2.0dev r12699" GEOS="3.5.0dev-CAPI-1.9.0 r3989" SFCGAL="1.0.4" PROJ="Rel. 4.8.0, 6 March 2012" -GDAL="GDAL 1.11.0, released 2014/04/16" LIBXML="2.7.8" LIBJSON="0.12" RASTER -(1 row) +POSTGIS="3.0.0dev r17211" [EXTENSION] PGSQL="110" GEOS="3.8.0dev-CAPI-1.11.0 df24b6bb" SFCGAL="1.3.6" PROJ="Rel. 5.2.0, September 15th, 2018" +GDAL="GDAL 2.3.2, released 2018/09/21" LIBXML="2.9.9" LIBJSON="0.13.1" LIBPROTOBUF="1.3.1" WAGYU="0.4.3 (Internal)" TOPOLOGY RASTER +(1 row) + @@ -465,6 +466,7 @@ GDAL="GDAL 1.11.0, released 2014/04/16" LIBXML="2.7.8" LIBJSON="0.12" RASTER , , , + , @@ -732,6 +734,52 @@ postgis_liblwgeom_version + + + PostGIS_Wagyu_Version + + Returns the version number of the internal Wagyu library. + + + + + + text PostGIS_Wagyu_Version + + + + + + + + Description + + Returns the version number of the internal Wagyu library, or + NULL if Wagyu support is not enabled. + + + + Examples + + SELECT PostGIS_Wagyu_Version(); + postgis_wagyu_version +----------------------- + 0.4.3 (Internal) +(1 row) + + + + See Also + + , , , , , + + + PostGIS_Scripts_Build_Date diff --git a/doc/reference_output.xml b/doc/reference_output.xml index 28f53b1e5..49bf03705 100644 --- a/doc/reference_output.xml +++ b/doc/reference_output.xml @@ -1447,6 +1447,10 @@ SELECT ST_GeoHash(ST_SetSRID(ST_MakePoint(-126,48),4326),5); clip_geom is a boolean to control if geometries should be clipped or encoded as is. If NULL it will default to true. Availability: 2.4.0 + + + From 3.0, Wagyu can be chosen at configure time to clip and validate MVT polygons. This library is faster and produces more correct results than the GEOS default, but it might drop small polygons. + @@ -1457,11 +1461,19 @@ SELECT ST_GeoHash(ST_SetSRID(ST_MakePoint(-126,48),4326),5); 4096, 0, false)); st_astext -------------------------------------------------------------------- - MULTIPOLYGON(((5 4096,10 4096,10 4091,5 4096)),((5 4096,0 4096,0 4101,5 4096))) + MULTIPOLYGON(((5 4096,10 4091,10 4096,5 4096)),((5 4096,0 4101,0 4096,5 4096))) ]]> + + + See Also + + , + + + @@ -1543,7 +1555,7 @@ SELECT ST_GeoHash(ST_SetSRID(ST_MakePoint(-126,48),4326),5); ST_MakeBox2D(ST_Point(0, 0), ST_Point(4096, 4096)), 4096, 0, false) AS geom) AS q; st_asmvt -------------------------------------------------------------------- - \x1a320a0474657374121d12020000180322150946ec3f1a14453b0a09280f091413121e09091e0f1a026331220228012880207802 + \x1a320a0474657374121d1202000018032215095aa63f1a134631130a270f09280a121d0a14140f1a026331220228012880207802 ]]> diff --git a/macros/ax_cxx_compile_stdcxx.m4 b/macros/ax_cxx_compile_stdcxx.m4 new file mode 100644 index 000000000..9e9eaedaa --- /dev/null +++ b/macros/ax_cxx_compile_stdcxx.m4 @@ -0,0 +1,948 @@ +# =========================================================================== +# https://www.gnu.org/software/autoconf-archive/ax_cxx_compile_stdcxx.html +# =========================================================================== +# +# SYNOPSIS +# +# AX_CXX_COMPILE_STDCXX(VERSION, [ext|noext], [mandatory|optional]) +# +# DESCRIPTION +# +# Check for baseline language coverage in the compiler for the specified +# version of the C++ standard. If necessary, add switches to CXX and +# CXXCPP to enable support. VERSION may be '11' (for the C++11 standard) +# or '14' (for the C++14 standard). +# +# The second argument, if specified, indicates whether you insist on an +# extended mode (e.g. -std=gnu++11) or a strict conformance mode (e.g. +# -std=c++11). If neither is specified, you get whatever works, with +# preference for an extended mode. +# +# The third argument, if specified 'mandatory' or if left unspecified, +# indicates that baseline support for the specified C++ standard is +# required and that the macro should error out if no mode with that +# support is found. If specified 'optional', then configuration proceeds +# regardless, after defining HAVE_CXX${VERSION} if and only if a +# supporting mode is found. +# +# LICENSE +# +# Copyright (c) 2008 Benjamin Kosnik +# Copyright (c) 2012 Zack Weinberg +# Copyright (c) 2013 Roy Stogner +# Copyright (c) 2014, 2015 Google Inc.; contributed by Alexey Sokolov +# Copyright (c) 2015 Paul Norman +# Copyright (c) 2015 Moritz Klammler +# Copyright (c) 2016, 2018 Krzesimir Nowak +# +# Copying and distribution of this file, with or without modification, are +# permitted in any medium without royalty provided the copyright notice +# and this notice are preserved. This file is offered as-is, without any +# warranty. + +#serial 10 + +dnl This macro is based on the code from the AX_CXX_COMPILE_STDCXX_11 macro +dnl (serial version number 13). + +AC_DEFUN([AX_CXX_COMPILE_STDCXX], [dnl + m4_if([$1], [11], [ax_cxx_compile_alternatives="11 0x"], + [$1], [14], [ax_cxx_compile_alternatives="14 1y"], + [$1], [17], [ax_cxx_compile_alternatives="17 1z"], + [m4_fatal([invalid first argument `$1' to AX_CXX_COMPILE_STDCXX])])dnl + m4_if([$2], [], [], + [$2], [ext], [], + [$2], [noext], [], + [m4_fatal([invalid second argument `$2' to AX_CXX_COMPILE_STDCXX])])dnl + m4_if([$3], [], [ax_cxx_compile_cxx$1_required=true], + [$3], [mandatory], [ax_cxx_compile_cxx$1_required=true], + [$3], [optional], [ax_cxx_compile_cxx$1_required=false], + [m4_fatal([invalid third argument `$3' to AX_CXX_COMPILE_STDCXX])]) + AC_LANG_PUSH([C++])dnl + ac_success=no + + m4_if([$2], [noext], [], [dnl + if test x$ac_success = xno; then + for alternative in ${ax_cxx_compile_alternatives}; do + switch="-std=gnu++${alternative}" + cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_$switch]) + AC_CACHE_CHECK(whether $CXX supports C++$1 features with $switch, + $cachevar, + [ac_save_CXX="$CXX" + CXX="$CXX $switch" + AC_COMPILE_IFELSE([AC_LANG_SOURCE([_AX_CXX_COMPILE_STDCXX_testbody_$1])], + [eval $cachevar=yes], + [eval $cachevar=no]) + CXX="$ac_save_CXX"]) + if eval test x\$$cachevar = xyes; then + CXX="$CXX $switch" + if test -n "$CXXCPP" ; then + CXXCPP="$CXXCPP $switch" + fi + ac_success=yes + break + fi + done + fi]) + + m4_if([$2], [ext], [], [dnl + if test x$ac_success = xno; then + dnl HP's aCC needs +std=c++11 according to: + dnl http://h21007.www2.hp.com/portal/download/files/unprot/aCxx/PDF_Release_Notes/769149-001.pdf + dnl Cray's crayCC needs "-h std=c++11" + for alternative in ${ax_cxx_compile_alternatives}; do + for switch in -std=c++${alternative} +std=c++${alternative} "-h std=c++${alternative}"; do + cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_$switch]) + AC_CACHE_CHECK(whether $CXX supports C++$1 features with $switch, + $cachevar, + [ac_save_CXX="$CXX" + CXX="$CXX $switch" + AC_COMPILE_IFELSE([AC_LANG_SOURCE([_AX_CXX_COMPILE_STDCXX_testbody_$1])], + [eval $cachevar=yes], + [eval $cachevar=no]) + CXX="$ac_save_CXX"]) + if eval test x\$$cachevar = xyes; then + CXX="$CXX $switch" + if test -n "$CXXCPP" ; then + CXXCPP="$CXXCPP $switch" + fi + ac_success=yes + break + fi + done + if test x$ac_success = xyes; then + break + fi + done + fi]) + AC_LANG_POP([C++]) + if test x$ax_cxx_compile_cxx$1_required = xtrue; then + if test x$ac_success = xno; then + AC_MSG_ERROR([*** A compiler with support for C++$1 language features is required.]) + fi + fi + if test x$ac_success = xno; then + HAVE_CXX$1=0 + AC_MSG_NOTICE([No compiler with C++$1 support was found]) + else + HAVE_CXX$1=1 + AC_DEFINE(HAVE_CXX$1,1, + [define if the compiler supports basic C++$1 syntax]) + fi + AC_SUBST(HAVE_CXX$1) +]) + + +dnl Test body for checking C++11 support + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_11], + _AX_CXX_COMPILE_STDCXX_testbody_new_in_11 +) + + +dnl Test body for checking C++14 support + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_14], + _AX_CXX_COMPILE_STDCXX_testbody_new_in_11 + _AX_CXX_COMPILE_STDCXX_testbody_new_in_14 +) + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_17], + _AX_CXX_COMPILE_STDCXX_testbody_new_in_11 + _AX_CXX_COMPILE_STDCXX_testbody_new_in_14 + _AX_CXX_COMPILE_STDCXX_testbody_new_in_17 +) + +dnl Tests for new features in C++11 + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_11], [[ + +// If the compiler admits that it is not ready for C++11, why torture it? +// Hopefully, this will speed up the test. + +#ifndef __cplusplus + +#error "This is not a C++ compiler" + +#elif __cplusplus < 201103L + +#error "This is not a C++11 compiler" + +#else + +namespace cxx11 +{ + + namespace test_static_assert + { + + template + struct check + { + static_assert(sizeof(int) <= sizeof(T), "not big enough"); + }; + + } + + namespace test_final_override + { + + struct Base + { + virtual void f() {} + }; + + struct Derived : public Base + { + virtual void f() override {} + }; + + } + + namespace test_double_right_angle_brackets + { + + template < typename T > + struct check {}; + + typedef check single_type; + typedef check> double_type; + typedef check>> triple_type; + typedef check>>> quadruple_type; + + } + + namespace test_decltype + { + + int + f() + { + int a = 1; + decltype(a) b = 2; + return a + b; + } + + } + + namespace test_type_deduction + { + + template < typename T1, typename T2 > + struct is_same + { + static const bool value = false; + }; + + template < typename T > + struct is_same + { + static const bool value = true; + }; + + template < typename T1, typename T2 > + auto + add(T1 a1, T2 a2) -> decltype(a1 + a2) + { + return a1 + a2; + } + + int + test(const int c, volatile int v) + { + static_assert(is_same::value == true, ""); + static_assert(is_same::value == false, ""); + static_assert(is_same::value == false, ""); + auto ac = c; + auto av = v; + auto sumi = ac + av + 'x'; + auto sumf = ac + av + 1.0; + static_assert(is_same::value == true, ""); + static_assert(is_same::value == true, ""); + static_assert(is_same::value == true, ""); + static_assert(is_same::value == false, ""); + static_assert(is_same::value == true, ""); + return (sumf > 0.0) ? sumi : add(c, v); + } + + } + + namespace test_noexcept + { + + int f() { return 0; } + int g() noexcept { return 0; } + + static_assert(noexcept(f()) == false, ""); + static_assert(noexcept(g()) == true, ""); + + } + + namespace test_constexpr + { + + template < typename CharT > + unsigned long constexpr + strlen_c_r(const CharT *const s, const unsigned long acc) noexcept + { + return *s ? strlen_c_r(s + 1, acc + 1) : acc; + } + + template < typename CharT > + unsigned long constexpr + strlen_c(const CharT *const s) noexcept + { + return strlen_c_r(s, 0UL); + } + + static_assert(strlen_c("") == 0UL, ""); + static_assert(strlen_c("1") == 1UL, ""); + static_assert(strlen_c("example") == 7UL, ""); + static_assert(strlen_c("another\0example") == 7UL, ""); + + } + + namespace test_rvalue_references + { + + template < int N > + struct answer + { + static constexpr int value = N; + }; + + answer<1> f(int&) { return answer<1>(); } + answer<2> f(const int&) { return answer<2>(); } + answer<3> f(int&&) { return answer<3>(); } + + void + test() + { + int i = 0; + const int c = 0; + static_assert(decltype(f(i))::value == 1, ""); + static_assert(decltype(f(c))::value == 2, ""); + static_assert(decltype(f(0))::value == 3, ""); + } + + } + + namespace test_uniform_initialization + { + + struct test + { + static const int zero {}; + static const int one {1}; + }; + + static_assert(test::zero == 0, ""); + static_assert(test::one == 1, ""); + + } + + namespace test_lambdas + { + + void + test1() + { + auto lambda1 = [](){}; + auto lambda2 = lambda1; + lambda1(); + lambda2(); + } + + int + test2() + { + auto a = [](int i, int j){ return i + j; }(1, 2); + auto b = []() -> int { return '0'; }(); + auto c = [=](){ return a + b; }(); + auto d = [&](){ return c; }(); + auto e = [a, &b](int x) mutable { + const auto identity = [](int y){ return y; }; + for (auto i = 0; i < a; ++i) + a += b--; + return x + identity(a + b); + }(0); + return a + b + c + d + e; + } + + int + test3() + { + const auto nullary = [](){ return 0; }; + const auto unary = [](int x){ return x; }; + using nullary_t = decltype(nullary); + using unary_t = decltype(unary); + const auto higher1st = [](nullary_t f){ return f(); }; + const auto higher2nd = [unary](nullary_t f1){ + return [unary, f1](unary_t f2){ return f2(unary(f1())); }; + }; + return higher1st(nullary) + higher2nd(nullary)(unary); + } + + } + + namespace test_variadic_templates + { + + template + struct sum; + + template + struct sum + { + static constexpr auto value = N0 + sum::value; + }; + + template <> + struct sum<> + { + static constexpr auto value = 0; + }; + + static_assert(sum<>::value == 0, ""); + static_assert(sum<1>::value == 1, ""); + static_assert(sum<23>::value == 23, ""); + static_assert(sum<1, 2>::value == 3, ""); + static_assert(sum<5, 5, 11>::value == 21, ""); + static_assert(sum<2, 3, 5, 7, 11, 13>::value == 41, ""); + + } + + // http://stackoverflow.com/questions/13728184/template-aliases-and-sfinae + // Clang 3.1 fails with headers of libstd++ 4.8.3 when using std::function + // because of this. + namespace test_template_alias_sfinae + { + + struct foo {}; + + template + using member = typename T::member_type; + + template + void func(...) {} + + template + void func(member*) {} + + void test(); + + void test() { func(0); } + + } + +} // namespace cxx11 + +#endif // __cplusplus >= 201103L + +]]) + + +dnl Tests for new features in C++14 + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_14], [[ + +// If the compiler admits that it is not ready for C++14, why torture it? +// Hopefully, this will speed up the test. + +#ifndef __cplusplus + +#error "This is not a C++ compiler" + +#elif __cplusplus < 201402L + +#error "This is not a C++14 compiler" + +#else + +namespace cxx14 +{ + + namespace test_polymorphic_lambdas + { + + int + test() + { + const auto lambda = [](auto&&... args){ + const auto istiny = [](auto x){ + return (sizeof(x) == 1UL) ? 1 : 0; + }; + const int aretiny[] = { istiny(args)... }; + return aretiny[0]; + }; + return lambda(1, 1L, 1.0f, '1'); + } + + } + + namespace test_binary_literals + { + + constexpr auto ivii = 0b0000000000101010; + static_assert(ivii == 42, "wrong value"); + + } + + namespace test_generalized_constexpr + { + + template < typename CharT > + constexpr unsigned long + strlen_c(const CharT *const s) noexcept + { + auto length = 0UL; + for (auto p = s; *p; ++p) + ++length; + return length; + } + + static_assert(strlen_c("") == 0UL, ""); + static_assert(strlen_c("x") == 1UL, ""); + static_assert(strlen_c("test") == 4UL, ""); + static_assert(strlen_c("another\0test") == 7UL, ""); + + } + + namespace test_lambda_init_capture + { + + int + test() + { + auto x = 0; + const auto lambda1 = [a = x](int b){ return a + b; }; + const auto lambda2 = [a = lambda1(x)](){ return a; }; + return lambda2(); + } + + } + + namespace test_digit_separators + { + + constexpr auto ten_million = 100'000'000; + static_assert(ten_million == 100000000, ""); + + } + + namespace test_return_type_deduction + { + + auto f(int& x) { return x; } + decltype(auto) g(int& x) { return x; } + + template < typename T1, typename T2 > + struct is_same + { + static constexpr auto value = false; + }; + + template < typename T > + struct is_same + { + static constexpr auto value = true; + }; + + int + test() + { + auto x = 0; + static_assert(is_same::value, ""); + static_assert(is_same::value, ""); + return x; + } + + } + +} // namespace cxx14 + +#endif // __cplusplus >= 201402L + +]]) + + +dnl Tests for new features in C++17 + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_17], [[ + +// If the compiler admits that it is not ready for C++17, why torture it? +// Hopefully, this will speed up the test. + +#ifndef __cplusplus + +#error "This is not a C++ compiler" + +#elif __cplusplus < 201703L + +#error "This is not a C++17 compiler" + +#else + +#include +#include +#include + +namespace cxx17 +{ + + namespace test_constexpr_lambdas + { + + constexpr int foo = [](){return 42;}(); + + } + + namespace test::nested_namespace::definitions + { + + } + + namespace test_fold_expression + { + + template + int multiply(Args... args) + { + return (args * ... * 1); + } + + template + bool all(Args... args) + { + return (args && ...); + } + + } + + namespace test_extended_static_assert + { + + static_assert (true); + + } + + namespace test_auto_brace_init_list + { + + auto foo = {5}; + auto bar {5}; + + static_assert(std::is_same, decltype(foo)>::value); + static_assert(std::is_same::value); + } + + namespace test_typename_in_template_template_parameter + { + + template typename X> struct D; + + } + + namespace test_fallthrough_nodiscard_maybe_unused_attributes + { + + int f1() + { + return 42; + } + + [[nodiscard]] int f2() + { + [[maybe_unused]] auto unused = f1(); + + switch (f1()) + { + case 17: + f1(); + [[fallthrough]]; + case 42: + f1(); + } + return f1(); + } + + } + + namespace test_extended_aggregate_initialization + { + + struct base1 + { + int b1, b2 = 42; + }; + + struct base2 + { + base2() { + b3 = 42; + } + int b3; + }; + + struct derived : base1, base2 + { + int d; + }; + + derived d1 {{1, 2}, {}, 4}; // full initialization + derived d2 {{}, {}, 4}; // value-initialized bases + + } + + namespace test_general_range_based_for_loop + { + + struct iter + { + int i; + + int& operator* () + { + return i; + } + + const int& operator* () const + { + return i; + } + + iter& operator++() + { + ++i; + return *this; + } + }; + + struct sentinel + { + int i; + }; + + bool operator== (const iter& i, const sentinel& s) + { + return i.i == s.i; + } + + bool operator!= (const iter& i, const sentinel& s) + { + return !(i == s); + } + + struct range + { + iter begin() const + { + return {0}; + } + + sentinel end() const + { + return {5}; + } + }; + + void f() + { + range r {}; + + for (auto i : r) + { + [[maybe_unused]] auto v = i; + } + } + + } + + namespace test_lambda_capture_asterisk_this_by_value + { + + struct t + { + int i; + int foo() + { + return [*this]() + { + return i; + }(); + } + }; + + } + + namespace test_enum_class_construction + { + + enum class byte : unsigned char + {}; + + byte foo {42}; + + } + + namespace test_constexpr_if + { + + template + int f () + { + if constexpr(cond) + { + return 13; + } + else + { + return 42; + } + } + + } + + namespace test_selection_statement_with_initializer + { + + int f() + { + return 13; + } + + int f2() + { + if (auto i = f(); i > 0) + { + return 3; + } + + switch (auto i = f(); i + 4) + { + case 17: + return 2; + + default: + return 1; + } + } + + } + + namespace test_template_argument_deduction_for_class_templates + { + + template + struct pair + { + pair (T1 p1, T2 p2) + : m1 {p1}, + m2 {p2} + {} + + T1 m1; + T2 m2; + }; + + void f() + { + [[maybe_unused]] auto p = pair{13, 42u}; + } + + } + + namespace test_non_type_auto_template_parameters + { + + template + struct B + {}; + + B<5> b1; + B<'a'> b2; + + } + + namespace test_structured_bindings + { + + int arr[2] = { 1, 2 }; + std::pair pr = { 1, 2 }; + + auto f1() -> int(&)[2] + { + return arr; + } + + auto f2() -> std::pair& + { + return pr; + } + + struct S + { + int x1 : 2; + volatile double y1; + }; + + S f3() + { + return {}; + } + + auto [ x1, y1 ] = f1(); + auto& [ xr1, yr1 ] = f1(); + auto [ x2, y2 ] = f2(); + auto& [ xr2, yr2 ] = f2(); + const auto [ x3, y3 ] = f3(); + + } + + namespace test_exception_spec_type_system + { + + struct Good {}; + struct Bad {}; + + void g1() noexcept; + void g2(); + + template + Bad + f(T*, T*); + + template + Good + f(T1*, T2*); + + static_assert (std::is_same_v); + + } + + namespace test_inline_variables + { + + template void f(T) + {} + + template inline T g(T) + { + return T{}; + } + + template<> inline void f<>(int) + {} + + template<> int g<>(int) + { + return 5; + } + + } + +} // namespace cxx17 + +#endif // __cplusplus < 201703L + +]]) diff --git a/postgis/Makefile.in b/postgis/Makefile.in index 7e975f3b3..3fce06d32 100644 --- a/postgis/Makefile.in +++ b/postgis/Makefile.in @@ -60,6 +60,13 @@ ifeq (@HAVE_PROTOBUF@,yes) PROTOBUF_OBJ=vector_tile.pb-c.o geobuf.pb-c.o endif +WAGYU_LIBPATH = ../deps/wagyu/@WAGYU_LIB@ +ifeq (@HAVE_WAGYU@,yes) +WAYGU_INCLUDE += -I../deps/wagyu +WAYGU_LIB = $(WAGYU_LIBPATH) +WAGYU_DEP = $(WAGYU_LIBPATH) +endif + # SQL preprocessor SQLPP = @SQLPP@ @@ -129,8 +136,8 @@ OBJS=$(PG_OBJS) # to an existing liblwgeom.so in the PostgreSQL $libdir supplied by an # older version of PostGIS, rather than with the static liblwgeom.a # supplied with newer versions of PostGIS -PG_CPPFLAGS += -I../liblwgeom @CFLAGS@ -I../libpgcommon @CPPFLAGS@ -fPIC -SHLIB_LINK_F = ../libpgcommon/libpgcommon.a ../liblwgeom/.libs/liblwgeom.a @SHLIB_LINK@ +PG_CPPFLAGS += -I../liblwgeom @CFLAGS@ -I../libpgcommon $(WAYGU_INCLUDE) @CPPFLAGS@ -fPIC +SHLIB_LINK_F = $(WAYGU_LIB) ../libpgcommon/libpgcommon.a ../liblwgeom/.libs/liblwgeom.a @SHLIB_LINK@ # Add SFCGAL Flags if defined ifeq (@SFCGAL@,sfcgal) @@ -183,7 +190,7 @@ endif # Also they are all dependent on postgis_config.h # and thus postgis_svn_revision.h # -$(PG_OBJS): ../liblwgeom/.libs/liblwgeom.a ../libpgcommon/libpgcommon.a ../postgis_config.h ../postgis_svn_revision.h +$(PG_OBJS): ../liblwgeom/.libs/liblwgeom.a ../libpgcommon/libpgcommon.a ../postgis_config.h ../postgis_svn_revision.h $(WAGYU_DEP) ../postgis_svn_revision.h: $(MAKE) -C .. postgis_svn_revision.h @@ -194,6 +201,9 @@ $(PG_OBJS): ../liblwgeom/.libs/liblwgeom.a ../libpgcommon/libpgcommon.a ../postg ../libpgcommon/libpgcommon.a: $(MAKE) -C ../libpgcommon libpgcommon.a +$(WAGYU_LIBPATH): + $(MAKE) -C ../deps/wagyu all + # Get protobuf-c compiler command PROTOCC=@PROTOCC@ diff --git a/postgis/mvt.c b/postgis/mvt.c index 2190e28f0..5c4be1eb1 100644 --- a/postgis/mvt.c +++ b/postgis/mvt.c @@ -553,7 +553,6 @@ static void parse_datum_as_string(mvt_agg_context *ctx, Oid typoid, add_value_as_string(ctx, value, tags, k); } -#if POSTGIS_PGSQL_VERSION >= 94 static uint32_t *parse_jsonb(mvt_agg_context *ctx, Jsonb *jb, uint32_t *tags) { @@ -631,7 +630,6 @@ static uint32_t *parse_jsonb(mvt_agg_context *ctx, Jsonb *jb, return tags; } -#endif /** * Sets the feature id. Ignores Nulls and negative values @@ -914,6 +912,48 @@ mvt_clip_and_validate_geos(LWGEOM *lwgeom, uint8_t basic_type, uint32_t extent, return ng; } +#ifdef HAVE_WAGYU + +#include "lwgeom_wagyu.h" + +static LWGEOM * +mvt_clip_and_validate(LWGEOM *lwgeom, uint8_t basic_type, uint32_t extent, uint32_t buffer, bool clip_geom) +{ + GBOX clip_box = {0}; + LWGEOM *clipped_lwgeom; + + /* Wagyu only supports polygons. Default to geos for other types */ + lwgeom = lwgeom_to_basic_type(lwgeom, POLYGONTYPE); + if (lwgeom->type != POLYGONTYPE && lwgeom->type != MULTIPOLYGONTYPE) + { + return mvt_clip_and_validate_geos(lwgeom, basic_type, extent, buffer, clip_geom); + } + + if (!clip_geom) + { + /* With clipping disabled, we request a clip with the geometry bbox to force validation */ + lwgeom_calculate_gbox(lwgeom, &clip_box); + } + else + { + clip_box.xmax = clip_box.ymax = (double)extent + (double)buffer; + clip_box.xmin = clip_box.ymin = -(double)buffer; + } + + clipped_lwgeom = lwgeom_wagyu_clip_by_box(lwgeom, &clip_box); + + return clipped_lwgeom; +} + +#else /* ! HAVE_WAGYU */ + +static LWGEOM * +mvt_clip_and_validate(LWGEOM *lwgeom, uint8_t basic_type, uint32_t extent, uint32_t buffer, bool clip_geom) +{ + return mvt_clip_and_validate_geos(lwgeom, basic_type, extent, buffer, clip_geom); +} +#endif + /** * Transform a geometry into vector tile coordinate space. * @@ -963,13 +1003,13 @@ LWGEOM *mvt_geom(LWGEOM *lwgeom, const GBOX *gbox, uint32_t extent, uint32_t buf affine.yoff = -gbox->ymax * fy; lwgeom_affine(lwgeom, &affine); - /* snap to integer precision, removing duplicate points */ + /* Snap to integer precision, removing duplicate points */ lwgeom_grid_in_place(lwgeom, &grid); if (lwgeom == NULL || lwgeom_is_empty(lwgeom)) return NULL; - lwgeom = mvt_clip_and_validate_geos(lwgeom, basic_type, extent, buffer, clip_geom); + lwgeom = mvt_clip_and_validate(lwgeom, basic_type, extent, buffer, clip_geom); if (lwgeom == NULL || lwgeom_is_empty(lwgeom)) return NULL; diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in index d59be49f0..e015e324a 100644 --- a/postgis/postgis.sql.in +++ b/postgis/postgis.sql.in @@ -2808,6 +2808,10 @@ CREATE OR REPLACE FUNCTION postgis_proj_version() RETURNS text AS 'MODULE_PATHNAME' LANGUAGE 'c' IMMUTABLE; +CREATE OR REPLACE FUNCTION postgis_wagyu_version() RETURNS text + AS 'MODULE_PATHNAME' + LANGUAGE 'c' IMMUTABLE; + -- -- IMPORTANT: -- Starting at 1.1.0 this function is used by postgis_proc_upgrade.pl @@ -2963,6 +2967,7 @@ DECLARE topo_scr_ver text := NULL; json_lib_ver text; protobuf_lib_ver text; + wagyu_lib_ver text; sfcgal_lib_ver text; sfcgal_scr_ver text; pgsql_scr_ver text; @@ -2974,6 +2979,7 @@ BEGIN SELECT @extschema@.postgis_geos_version() INTO geosver; SELECT @extschema@.postgis_libjson_version() INTO json_lib_ver; SELECT @extschema@.postgis_libprotobuf_version() INTO protobuf_lib_ver; + SELECT @extschema@.postgis_wagyu_version() INTO wagyu_lib_ver; SELECT @extschema@._postgis_scripts_pgsql_version() INTO pgsql_scr_ver; SELECT @extschema@._postgis_pgsql_version() INTO pgsql_ver; BEGIN @@ -3083,6 +3089,10 @@ BEGIN fullver = fullver || ' LIBPROTOBUF="' || protobuf_lib_ver || '"'; END IF; + IF wagyu_lib_ver IS NOT NULL THEN + fullver = fullver || ' WAGYU="' || wagyu_lib_ver || '"'; + END IF; + IF dbproc != relproc THEN fullver = fullver || ' (core procs from "' || dbproc || '" need upgrade)'; END IF; diff --git a/postgis/postgis_libprotobuf.c b/postgis/postgis_libprotobuf.c index 4f3a2ae5b..99e208fb1 100644 --- a/postgis/postgis_libprotobuf.c +++ b/postgis/postgis_libprotobuf.c @@ -7,6 +7,10 @@ #include #endif +#ifdef HAVE_WAGYU +#include "lwgeom_wagyu.h" +#endif + PG_FUNCTION_INFO_V1(postgis_libprotobuf_version); Datum postgis_libprotobuf_version(PG_FUNCTION_ARGS) { @@ -17,4 +21,16 @@ Datum postgis_libprotobuf_version(PG_FUNCTION_ARGS) text *result = cstring_to_text(ver); PG_RETURN_POINTER(result); #endif +} + +PG_FUNCTION_INFO_V1(postgis_wagyu_version); +Datum postgis_wagyu_version(PG_FUNCTION_ARGS) +{ +#ifndef HAVE_WAGYU + PG_RETURN_NULL(); +#else /* HAVE_WAGYU */ + const char *ver = libwagyu_version(); + text *result = cstring_to_text(ver); + PG_RETURN_POINTER(result); +#endif } \ No newline at end of file diff --git a/postgis/postgis_module.c b/postgis/postgis_module.c index 6631c7617..497bf7909 100644 --- a/postgis/postgis_module.c +++ b/postgis/postgis_module.c @@ -36,6 +36,10 @@ #include "geos_c.h" #include "lwgeom_backend_api.h" +#ifdef HAVE_WAGYU +#include "lwgeom_wagyu.h" +#endif + /* * This is required for builds against pgsql */ @@ -97,6 +101,10 @@ handleInterrupt(int sig) GEOS_interruptRequest(); +#ifdef HAVE_WAGYU + lwgeom_wagyu_interruptRequest(); +#endif + /* request interruption of liblwgeom as well */ lwgeom_request_interrupt(); diff --git a/postgis_config.h.in b/postgis_config.h.in index 0b14758fd..15ed772cb 100644 --- a/postgis_config.h.in +++ b/postgis_config.h.in @@ -110,6 +110,9 @@ /* Define to 1 if you have the header file. */ #undef HAVE_UNISTD_H +/* Define to 1 if wagyu is being built */ +#undef HAVE_WAGYU + /* Define to the sub-directory in which libtool stores uninstalled libraries. */ #undef LT_OBJDIR diff --git a/regress/core/mvt.sql b/regress/core/mvt.sql index edbe1548c..b71ea39e1 100644 --- a/regress/core/mvt.sql +++ b/regress/core/mvt.sql @@ -31,14 +31,20 @@ select 'PG6', ST_AsText(ST_AsMVTGeom( ST_MakeBox2D(ST_Point(626172.135625, 6261721.35625), ST_Point(1252344.27125, 6887893.49188)), 4096, 0, false)); -SELECT 'PG7', ST_AsText(ST_AsMVTGeom( +WITH geometry AS +( +SELECT ST_AsText(ST_AsMVTGeom( ST_GeomFromText('POLYGON((-7792023.4539488 1411512.60791779,-7785283.40665468 1406282.69482469,-7783978.88137195 1404858.20373788,-7782986.89858399 1402324.91434802,-7779028.02672366 1397370.31802772, -7778652.06985644 1394387.75452545,-7779906.76953697 1393279.22658385,-7782212.33678782 1393293.14086794,-7784631.14401331 1394225.4151684,-7786257.27108231 1395867.40241344,-7783978.88137195 1395867.40241344, -7783978.88137195 1396646.68250521,-7787752.03959369 1398469.72134299,-7795443.30325373 1405280.43988858,-7797717.16326269 1406217.73286975,-7798831.44531677 1406904.48130551,-7799311.5830898 1408004.24038921, -7799085.10302919 1409159.72782477,-7798052.35381919 1411108.84582812,-7797789.63692662 1412213.40365339,-7798224.47868753 1414069.89725829,-7799003.5701851 1415694.42577482,-7799166.63587328 1416966.26267896, -7797789.63692662 1417736.81850415,-7793160.38395328 1412417.61222784,-7792023.4539488 1411512.60791779))'), ST_MakeBox2D(ST_Point(-20037508.34, -20037508.34), ST_Point(20037508.34, 20037508.34)), - 4096, 10, true)) as g; + 4096, 10, true)) as g +) +-- Wagyu might drop small polygons like this one bellow +-- See https://github.com/mapbox/wagyu/issues/94#issuecomment-452376133 +SELECT 'PG7', COALESCE(g, 'POLYGON((1252 1904,1253 1905,1253 1906,1251 1904,1252 1904))') FROM geometry; select 'PG8', ST_AsText(ST_AsMVTGeom( ST_GeomFromText('GEOMETRYCOLLECTION(MULTIPOLYGON (((0 0, 10 0, 10 5, 0 -5, 0 0))))'), @@ -49,7 +55,7 @@ select 'PG9', ST_Area(ST_AsMVTGeom( ST_MakeBox2D(ST_Point(0, 0), ST_Point(5, 5)), 4096, 0, true)); --- There shoulnd't be floating point values +-- There shouldn't be floating point values WITH geometry AS ( SELECT ST_AsMVTGeom( @@ -359,11 +365,17 @@ SELECT 'PG56', ST_AsText(ST_AsMVTGeom( ST_MakeBox2D(ST_Point(0, 0), ST_Point(100, 100)), 100, 0, true)); --- This clips in float -SELECT 'PG57', ST_AsText(ST_AsMVTGeom( +-- Different round behaviour between geos and wagyu +WITH geometry AS +( + SELECT ST_AsText(ST_AsMVTGeom( ST_GeomFromText('POLYGON((0 0, 0 99, 1 101, 100 100, 100 0, 0 0))'), ST_MakeBox2D(ST_Point(0, 0), ST_Point(100, 100)), - 100, 0, true)); + 100, 0, true)) as g +) +SELECT 'PG57', + g = 'POLYGON((100 0,100 100,0 100,0 1,1 0,100 0))' OR g = 'POLYGON((0 1,0 0,100 0,100 100,0 100,0 1))' +FROM geometry; -- Geometrycollection test SELECT 'PG58', ST_AsText(ST_AsMVTGeom( @@ -391,10 +403,18 @@ SELECT 'PG62', ST_AsText(ST_AsMVTGeom( ST_MakeBox2D(ST_Point(0, 0), ST_Point(100, 100)), 100, 0, true)); -SELECT 'PG63', ST_AsText(ST_AsMVTGeom( - ST_GeomFromText('GEOMETRYCOLLECTION(LINESTRING(10 10, 20 20), POLYGON((90 90, 110 90, 110 110, 90 110, 90 90)), LINESTRING(20 20, 15 15))'), +-- Same polygon, different starting point between Wagyu and GEOS backends +WITH geometry AS +( + SELECT ST_AsText(ST_AsMVTGeom( + ST_GeomFromText('GEOMETRYCOLLECTION(LINESTRING(10 10, 20 20), POLYGON((110 90, 110 110, 90 110, 90 90, 110 90)), LINESTRING(20 20, 15 15))'), ST_MakeBox2D(ST_Point(0, 0), ST_Point(100, 100)), - 100, 0, true)); + 100, 0, true)) as g +) +SELECT 'PG63', + ST_Area(g), + g = 'POLYGON((90 10,90 0,100 0,100 10,90 10))' OR g = 'POLYGON((90 0,100 0,100 10,90 10,90 0))' +FROM geometry; SELECT 'PG64', ST_AsText(ST_AsMVTGeom( ST_GeomFromText('GEOMETRYCOLLECTION(MULTIPOLYGON EMPTY, POINT(50 50))'), diff --git a/regress/core/mvt_expected b/regress/core/mvt_expected index 6aa2742e5..a4960cf45 100644 --- a/regress/core/mvt_expected +++ b/regress/core/mvt_expected @@ -65,13 +65,13 @@ PG53|600 PG54| PG55| PG56| -PG57|POLYGON((0 1,0 0,100 0,100 100,0 100,0 1)) +PG57|t PG58|POLYGON((0 100,0 90,10 90,10 100,0 100)) PG59|POLYGON((0 100,0 90,10 90,10 100,0 100)) PG60|POLYGON((0 100,0 90,10 90,10 100,0 100)) PG61| PG62|POLYGON((0 100,0 90,10 90,10 100,0 100)) -PG63|POLYGON((90 0,100 0,100 10,90 10,90 0)) +PG63|100|t PG64| TG1|GiEKBHRlc3QSDBICAAAYASIECTLePxoCYzEiAigBKIAgeAI= TG2|GiMKBHRlc3QSDhICAAAYASIGETLePwIBGgJjMSICKAEogCB4Ag== diff --git a/utils/postgis_restore.pl.in b/utils/postgis_restore.pl.in index c3bc57ccf..0127a7d10 100644 --- a/utils/postgis_restore.pl.in +++ b/utils/postgis_restore.pl.in @@ -648,6 +648,7 @@ COMMENT FUNCTION postgis_lib_build_date() COMMENT FUNCTION postgis_lib_version() COMMENT FUNCTION postgis_libxml_version() COMMENT FUNCTION postgis_proj_version() +COMMENT FUNCTION postgis_wagyu_version() COMMENT FUNCTION postgis_raster_lib_build_date() COMMENT FUNCTION postgis_raster_lib_version() COMMENT FUNCTION postgis_scripts_build_date() -- 2.40.0