diff --git a/autoremesher.pro b/autoremesher.pro index 3fb7e0fa..bcc6288c 100644 --- a/autoremesher.pro +++ b/autoremesher.pro @@ -114,7 +114,7 @@ HEADERS += src/AutoRemesher/meshseparator.h SOURCES += src/AutoRemesher/relativeheight.cpp HEADERS += src/AutoRemesher/relativeheight.h -SOURCES += src/tiny_obj_loader.cc +# SOURCES += src/tiny_obj_loader.cc HEADERS += src/tiny_obj_loader.h INCLUDEPATH += thirdparty/openvdb/openvdb-7.0.0 @@ -421,6 +421,9 @@ INCLUDEPATH += thirdparty/xatlas/source/xatlas HEADERS += thirdparty/xatlas/source/xatlas/xatlas.h SOURCES += thirdparty/xatlas/source/xatlas/xatlas.cpp +INCLUDEPATH += thirdparty/earcut/include +HEADERS += thirdparty/earcut/include/mapbox/earcut.hpp + INCLUDEPATH += thirdparty/tbb/include unix { LIBS += -Lthirdparty/tbb/build2 -ltbbmalloc_proxy_static -ltbbmalloc_static -ltbb_static diff --git a/src/main.cpp b/src/main.cpp index ac2cf5c0..87c1033f 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -4,11 +4,251 @@ #include #include #include +#include +#define TINYOBJLOADER_IMPLEMENTATION #include "tiny_obj_loader.h" #include #include #include +struct PolygonPoolItem +{ + std::size_t beginIdx = 0, count = 0; + std::vector original; +}; + +// Ripped from https://github.com/tinyobjloader/tinyobjloader/blob/release/tiny_obj_loader.h#L1482 +std::vector triangulate(std::vector const &v, std::vector const &polygonIndices) +{ + std::vector result; + + using namespace tinyobj; + + face_t face; + face.smoothing_group_id = 0; + for (std::size_t idx : polygonIndices) + { + face.vertex_indices.emplace_back(vertex_index_t(idx)); + } + size_t npolys = face.vertex_indices.size(); + + if (npolys == 4) + { + vertex_index_t i0 = face.vertex_indices[0]; + vertex_index_t i1 = face.vertex_indices[1]; + vertex_index_t i2 = face.vertex_indices[2]; + vertex_index_t i3 = face.vertex_indices[3]; + + size_t vi0 = size_t(i0.v_idx); + size_t vi1 = size_t(i1.v_idx); + size_t vi2 = size_t(i2.v_idx); + size_t vi3 = size_t(i3.v_idx); + + if (((3 * vi0 + 2) >= v.size()) || ((3 * vi1 + 2) >= v.size()) || + ((3 * vi2 + 2) >= v.size()) || ((3 * vi3 + 2) >= v.size())) + { + // Invalid triangle. + // FIXME(syoyo): Is it ok to simply skip this invalid triangle? + return std::vector(); + } + + real_t v0x = v[vi0 * 3 + 0]; + real_t v0y = v[vi0 * 3 + 1]; + real_t v0z = v[vi0 * 3 + 2]; + real_t v1x = v[vi1 * 3 + 0]; + real_t v1y = v[vi1 * 3 + 1]; + real_t v1z = v[vi1 * 3 + 2]; + real_t v2x = v[vi2 * 3 + 0]; + real_t v2y = v[vi2 * 3 + 1]; + real_t v2z = v[vi2 * 3 + 2]; + real_t v3x = v[vi3 * 3 + 0]; + real_t v3y = v[vi3 * 3 + 1]; + real_t v3z = v[vi3 * 3 + 2]; + + // There are two candidates to split the quad into two triangles. + // + // Choose the shortest edge. + // TODO: Is it better to determine the edge to split by calculating + // the area of each triangle? + // + // +---+ + // |\ | + // | \ | + // | \| + // +---+ + // + // +---+ + // | /| + // | / | + // |/ | + // +---+ + + real_t e02x = v2x - v0x; + real_t e02y = v2y - v0y; + real_t e02z = v2z - v0z; + real_t e13x = v3x - v1x; + real_t e13y = v3y - v1y; + real_t e13z = v3z - v1z; + + real_t sqr02 = e02x * e02x + e02y * e02y + e02z * e02z; + real_t sqr13 = e13x * e13x + e13y * e13y + e13z * e13z; + + index_t idx0, idx1, idx2, idx3; + + idx0.vertex_index = i0.v_idx; + idx0.normal_index = i0.vn_idx; + idx0.texcoord_index = i0.vt_idx; + idx1.vertex_index = i1.v_idx; + idx1.normal_index = i1.vn_idx; + idx1.texcoord_index = i1.vt_idx; + idx2.vertex_index = i2.v_idx; + idx2.normal_index = i2.vn_idx; + idx2.texcoord_index = i2.vt_idx; + idx3.vertex_index = i3.v_idx; + idx3.normal_index = i3.vn_idx; + idx3.texcoord_index = i3.vt_idx; + + if (sqr02 < sqr13) + { + // [0, 1, 2], [0, 2, 3] + result.emplace_back(idx0.vertex_index); + result.emplace_back(idx1.vertex_index); + result.emplace_back(idx2.vertex_index); + + result.emplace_back(idx0.vertex_index); + result.emplace_back(idx2.vertex_index); + result.emplace_back(idx3.vertex_index); + } + else + { + // [0, 1, 3], [1, 2, 3] + result.emplace_back(idx0.vertex_index); + result.emplace_back(idx1.vertex_index); + result.emplace_back(idx3.vertex_index); + + result.emplace_back(idx1.vertex_index); + result.emplace_back(idx2.vertex_index); + result.emplace_back(idx3.vertex_index); + } + } + else + { + vertex_index_t i0 = face.vertex_indices[0]; + vertex_index_t i0_2 = i0; + + // TMW change: Find the normal axis of the polygon using Newell's + // method + TinyObjPoint n; + for (size_t k = 0; k < npolys; ++k) + { + i0 = face.vertex_indices[k % npolys]; + size_t vi0 = size_t(i0.v_idx); + + size_t j = (k + 1) % npolys; + i0_2 = face.vertex_indices[j]; + size_t vi0_2 = size_t(i0_2.v_idx); + + real_t v0x = v[vi0 * 3 + 0]; + real_t v0y = v[vi0 * 3 + 1]; + real_t v0z = v[vi0 * 3 + 2]; + + real_t v0x_2 = v[vi0_2 * 3 + 0]; + real_t v0y_2 = v[vi0_2 * 3 + 1]; + real_t v0z_2 = v[vi0_2 * 3 + 2]; + + const TinyObjPoint point1(v0x, v0y, v0z); + const TinyObjPoint point2(v0x_2, v0y_2, v0z_2); + + TinyObjPoint a(point1.x - point2.x, point1.y - point2.y, + point1.z - point2.z); + TinyObjPoint b(point1.x + point2.x, point1.y + point2.y, + point1.z + point2.z); + + n.x += (a.y * b.z); + n.y += (a.z * b.x); + n.z += (a.x * b.y); + } + real_t length_n = GetLength(n); + // Check if zero length normal + if (length_n <= 0) + { + // Failure - ignore face + return std::vector(); + } + // Negative is to flip the normal to the correct direction + real_t inv_length = -real_t(1.0) / length_n; + n.x *= inv_length; + n.y *= inv_length; + n.z *= inv_length; + + TinyObjPoint axis_w, axis_v, axis_u; + axis_w = n; + TinyObjPoint a; + if (std::fabs(axis_w.x) > real_t(0.9999999)) + { + a = TinyObjPoint(0, 1, 0); + } + else + { + a = TinyObjPoint(1, 0, 0); + } + axis_v = Normalize(cross(axis_w, a)); + axis_u = cross(axis_w, axis_v); + using Point = std::array; + + // first polyline define the main polygon. + // following polylines define holes(not used in tinyobj). + std::vector> polygon; + + std::vector polyline; + + // TMW change: Find best normal and project v0x and v0y to those + // coordinates, instead of picking a plane aligned with an axis (which + // can flip polygons). + + // Fill polygon data(facevarying vertices). + for (size_t k = 0; k < npolys; k++) + { + i0 = face.vertex_indices[k]; + size_t vi0 = size_t(i0.v_idx); + + assert(((3 * vi0 + 2) < v.size())); + + real_t v0x = v[vi0 * 3 + 0]; + real_t v0y = v[vi0 * 3 + 1]; + real_t v0z = v[vi0 * 3 + 2]; + + TinyObjPoint polypoint(v0x, v0y, v0z); + TinyObjPoint loc = WorldToLocal(polypoint, axis_u, axis_v, axis_w); + + polyline.push_back({loc.x, loc.y}); + } + + polygon.push_back(polyline); + std::vector indices = mapbox::earcut(polygon); + // => result = 3 * faces, clockwise + + assert(indices.size() % 3 == 0); + + // Reconstruct vertex_index_t + for (size_t k = 0; k < indices.size() / 3; k++) + { + { + index_t idx0, idx1, idx2; + idx0.vertex_index = face.vertex_indices[indices[3 * k + 0]].v_idx; + idx1.vertex_index = face.vertex_indices[indices[3 * k + 1]].v_idx; + idx2.vertex_index = face.vertex_indices[indices[3 * k + 2]].v_idx; + + result.emplace_back(idx0.vertex_index); + result.emplace_back(idx1.vertex_index); + result.emplace_back(idx2.vertex_index); + } + } + } + + return result; +} + bool loadObj(std::string const &filename, std::vector &vertices, std::vector> &triangles) { tinyobj::attrib_t attributes; @@ -119,30 +359,59 @@ int main(int argc, char *argv[]) } } } - std::vector outIndices; - std::vector outFacesNumVertices; + std::vector polygonPool; + std::vector xatlasInputTris; + std::vector xatlasInputTrisKeep; for (std::vector>::const_iterator it = remeshedQuads.begin(); it != remeshedQuads.end(); ++it) { - std::uint8_t faceNumVertices = 0; - for (std::vector::const_iterator subIt = (*it).begin(); subIt != (*it).end(); ++subIt) + if (it->size() == 3) { - outIndices.emplace_back(remeshed2OutIdx[*subIt]); - ++faceNumVertices; + xatlasInputTris.emplace_back(remeshed2OutIdx[(*it)[0]]); + xatlasInputTrisKeep.emplace_back(true); + xatlasInputTris.emplace_back(remeshed2OutIdx[(*it)[1]]); + xatlasInputTrisKeep.emplace_back(true); + xatlasInputTris.emplace_back(remeshed2OutIdx[(*it)[2]]); + xatlasInputTrisKeep.emplace_back(true); + } + else if (it->size() > 3) + { + std::vector polygonIndices; + for (std::vector::const_iterator subIt = it->begin(); subIt != it->end(); ++subIt) + { + polygonIndices.emplace_back(remeshed2OutIdx[*subIt]); + } + + std::vector triangulatedIndices = triangulate(outVertices, polygonIndices); + xatlasInputTris.insert(xatlasInputTris.end(), triangulatedIndices.begin(), triangulatedIndices.end()); + for (std::size_t i = 0; i < triangulatedIndices.size(); ++i) + { + xatlasInputTrisKeep.emplace_back(true); + } + + // Only support quads in polygon pool for now + if (it->size() == 4) + { + PolygonPoolItem item; + item.beginIdx = xatlasInputTris.size() - triangulatedIndices.size(); + item.count = triangulatedIndices.size(); + item.original = polygonIndices; + polygonPool.emplace_back(item); + } + } + else + { + // Ignore degenerate faces } - outFacesNumVertices.emplace_back(faceNumVertices); } - // xatlas xatlas::Atlas *atlas = xatlas::Create(); xatlas::MeshDecl meshDecl; meshDecl.vertexCount = outVerticesCount; meshDecl.vertexPositionData = outVertices.data(); meshDecl.vertexPositionStride = sizeof(float) * 3; - meshDecl.indexCount = outIndices.size(); - meshDecl.indexData = outIndices.data(); + meshDecl.indexCount = xatlasInputTris.size(); + meshDecl.indexData = xatlasInputTris.data(); meshDecl.indexFormat = xatlas::IndexFormat::UInt32; - meshDecl.faceVertexCount = outFacesNumVertices.data(); - meshDecl.faceCount = outFacesNumVertices.size(); xatlas::AddMeshError error = xatlas::AddMesh(atlas, meshDecl, 1); if (error != xatlas::AddMeshError::Success) { @@ -157,7 +426,7 @@ int main(int argc, char *argv[]) std::cout << "Atlas resolution: " << atlas->width << "x" << atlas->height << "\n"; std::cout << "Total vertices after atlasing: " << atlas->meshes[0].vertexCount << "\n"; - // OBJ output + // OBJ Output: Vertices std::ofstream outFile(outFilename, std::ios::out); for (std::size_t i = 0; i < outVerticesCount; ++i) { @@ -168,18 +437,71 @@ int main(int argc, char *argv[]) const xatlas::Vertex &vertex = atlas->meshes[0].vertexArray[i]; outFile << "vt " << vertex.uv[0] / static_cast(atlas->width) << " " << vertex.uv[1] / static_cast(atlas->height) << "\n"; } - std::size_t totalFaceIndex = 0; - for (std::size_t i = 0; i < outFacesNumVertices.size(); ++i) + + // OBJ Output: Polygons + for (PolygonPoolItem const &item : polygonPool) { - outFile << "f"; - for (std::size_t j = 0; j < outFacesNumVertices[i]; ++j) + bool isValid = true; + + std::unordered_map old2New; + for (std::size_t i = item.beginIdx; i < item.beginIdx + item.count; ++i) + { + const std::size_t newIdx = atlas->meshes[0].indexArray[i]; + const xatlas::Vertex &xatlasVertex = atlas->meshes[0].vertexArray[newIdx]; + std::size_t oldIdx = xatlasVertex.xref; + if (old2New.find(oldIdx) == old2New.end()) + { + old2New[oldIdx] = newIdx; + } + else if (old2New[oldIdx] != newIdx) + { + isValid = false; + break; + } + } + + if (isValid) + { + outFile << "f"; + for (std::size_t i = 0; i < item.original.size(); ++i) + { + std::size_t oldIdx = item.original[i]; + const std::size_t newIdx = old2New[item.original[i]]; + outFile << " " << (oldIdx + 1) << "/" << (newIdx + 1); + } + outFile << "\n"; + + for (std::size_t i = item.beginIdx; i < item.beginIdx + item.count; ++i) + { + xatlasInputTrisKeep[i] = false; + } + } + } + + // OBJ Output: Triangles + for (std::size_t i = 0; i < xatlasInputTris.size(); i += 3) + { + bool wrote = false; + for (std::size_t j = 0; j < 3; ++j) + { + if (xatlasInputTrisKeep[i + j]) + { + if (!wrote) + { + outFile << "f"; + wrote = true; + } + + const std::size_t texIdx = atlas->meshes[0].indexArray[i + j]; + const xatlas::Vertex &xatlasVertex = atlas->meshes[0].vertexArray[texIdx]; + std::size_t posIdx = xatlasVertex.xref; + outFile << " " << (posIdx + 1) << "/" << (texIdx + 1); + } + } + if (wrote) { - const std::size_t index = atlas->meshes[0].indexArray[totalFaceIndex]; - const xatlas::Vertex &vertex = atlas->meshes[0].vertexArray[index]; - outFile << " " << (vertex.xref + 1) << "/" << (index + 1); - ++totalFaceIndex; + outFile << "\n"; } - outFile << "\n"; } std::cout << "AutoRemesher done!\n"; diff --git a/thirdparty/earcut/LICENSE b/thirdparty/earcut/LICENSE new file mode 100644 index 00000000..8bafb577 --- /dev/null +++ b/thirdparty/earcut/LICENSE @@ -0,0 +1,15 @@ +ISC License + +Copyright (c) 2015, 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. diff --git a/thirdparty/earcut/README.md b/thirdparty/earcut/README.md new file mode 100644 index 00000000..67f235b7 --- /dev/null +++ b/thirdparty/earcut/README.md @@ -0,0 +1,131 @@ +## Earcut + +A C++ port of [earcut.js](https://github.com/mapbox/earcut), a fast, [header-only](https://github.com/mapbox/earcut.hpp/blob/master/include/mapbox/earcut.hpp) polygon triangulation library. + +[![Travis](https://img.shields.io/travis/com/mapbox/earcut.hpp.svg)](https://travis-ci.com/github/mapbox/earcut.hpp) +[![AppVeyor](https://ci.appveyor.com/api/projects/status/a1ysrqd69mqn7coo/branch/master?svg=true)](https://ci.appveyor.com/project/Mapbox/earcut-hpp-8wm4o/branch/master) +[![Coverage](https://img.shields.io/coveralls/github/mapbox/earcut.hpp.svg)](https://coveralls.io/github/mapbox/earcut.hpp) +[![Coverity Scan](https://img.shields.io/coverity/scan/14000.svg)](https://scan.coverity.com/projects/14000) +[![Average time to resolve an issue](http://isitmaintained.com/badge/resolution/mapbox/earcut.hpp.svg)](http://isitmaintained.com/project/mapbox/earcut.hpp "Average time to resolve an issue") +[![Percentage of issues still open](http://isitmaintained.com/badge/open/mapbox/earcut.hpp.svg)](http://isitmaintained.com/project/mapbox/earcut.hpp "Percentage of issues still open") +[![Mourner](https://img.shields.io/badge/simply-awesome-brightgreen.svg)](https://github.com/mourner/projects) + +The library implements a modified ear slicing algorithm, optimized by [z-order curve](http://en.wikipedia.org/wiki/Z-order_curve) hashing and extended to handle holes, twisted polygons, degeneracies and self-intersections in a way that doesn't _guarantee_ correctness of triangulation, but attempts to always produce acceptable results for practical data like geographical shapes. + +It's based on ideas from [FIST: Fast Industrial-Strength Triangulation of Polygons](http://www.cosy.sbg.ac.at/~held/projects/triang/triang.html) by Martin Held and [Triangulation by Ear Clipping](http://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf) by David Eberly. + +## Usage + +```cpp +#include +``` +```cpp +// The number type to use for tessellation +using Coord = double; + +// The index type. Defaults to uint32_t, but you can also pass uint16_t if you know that your +// data won't have more than 65536 vertices. +using N = uint32_t; + +// Create array +using Point = std::array; +std::vector> polygon; + +// Fill polygon structure with actual data. Any winding order works. +// The first polyline defines the main polygon. +polygon.push_back({{100, 0}, {100, 100}, {0, 100}, {0, 0}}); +// Following polylines define holes. +polygon.push_back({{75, 25}, {75, 75}, {25, 75}, {25, 25}}); + +// Run tessellation +// Returns array of indices that refer to the vertices of the input polygon. +// e.g: the index 6 would refer to {25, 75} in this example. +// Three subsequent indices form a triangle. Output triangles are clockwise. +std::vector indices = mapbox::earcut(polygon); +``` + +Earcut can triangulate a simple, planar polygon of any winding order including holes. It will even return a robust, acceptable solution for non-simple poygons. Earcut works on a 2D plane. If you have three or more dimensions, you can project them onto a 2D surface before triangulation, or use a more suitable library for the task (e.g [CGAL](https://doc.cgal.org/latest/Triangulation_3/index.html)). + + +It is also possible to use your custom point type as input. There are default accessors defined for `std::tuple`, `std::pair`, and `std::array`. For a custom type (like Clipper's `IntPoint` type), do this: + +```cpp +// struct IntPoint { +// int64_t X, Y; +// }; + +namespace mapbox { +namespace util { + +template <> +struct nth<0, IntPoint> { + inline static auto get(const IntPoint &t) { + return t.X; + }; +}; +template <> +struct nth<1, IntPoint> { + inline static auto get(const IntPoint &t) { + return t.Y; + }; +}; + +} // namespace util +} // namespace mapbox +``` + +You can also use a custom container type for your polygon. Similar to std::vector, it has to meet the requirements of [Container](https://en.cppreference.com/w/cpp/named_req/Container), in particular `size()`, `empty()` and `operator[]`. + +

+ example triangulation +

+ +## Additional build instructions +In case you just want to use the earcut triangulation library; copy and include the header file [``](https://github.com/mapbox/earcut.hpp/blob/master/include/mapbox/earcut.hpp) in your project and follow the steps documented in the section [Usage](#usage). + +If you want to build the test, benchmark and visualization programs instead, follow these instructions: + +### Dependencies + +Before you continue, make sure to have the following tools and libraries installed: + * git ([Ubuntu](https://help.ubuntu.com/lts/serverguide/git.html)/[Windows/macOS](http://git-scm.com/downloads)) + * cmake 3.2+ ([Ubuntu](https://launchpad.net/~george-edison55/+archive/ubuntu/cmake-3.x)/[Windows/macOS](https://cmake.org/download/)) + * OpenGL SDK ([Ubuntu](http://packages.ubuntu.com/de/trusty/libgl1-mesa-dev)/[Windows](https://dev.windows.com/en-us/downloads/windows-10-sdk)/[macOS](https://developer.apple.com/opengl/)) + * Compiler such as [GCC 4.9+, Clang 3.4+](https://launchpad.net/~ubuntu-toolchain-r/+archive/ubuntu/test), [MSVC12+](https://www.visualstudio.com/) + +Note: On some operating systems such as Windows, manual steps are required to add cmake and [git](http://blog.countableset.ch/2012/06/07/adding-git-to-windows-7-path/) to your PATH environment variable. + +### Manual compilation + +```bash +git clone --recursive https://github.com/mapbox/earcut.hpp.git +cd earcut.hpp +mkdir build +cd build +cmake .. +make +# ./tests +# ./bench +# ./viz +``` + +### [Visual Studio](https://www.visualstudio.com/), [Eclipse](https://eclipse.org/), [XCode](https://developer.apple.com/xcode/), ... + +```batch +git clone --recursive https://github.com/mapbox/earcut.hpp.git +cd earcut.hpp +mkdir project +cd project +cmake .. -G "Visual Studio 14 2015" +::you can also generate projects for "Visual Studio 12 2013", "XCode", "Eclipse CDT4 - Unix Makefiles" +``` +After completion, open the generated project with your IDE. + + +### [CLion](https://www.jetbrains.com/clion/), [Visual Studio 2017+](https://www.visualstudio.com/) + +Import the project from https://github.com/mapbox/earcut.hpp.git and you should be good to go! + +## Status + +This is currently based on [earcut 2.2.4](https://github.com/mapbox/earcut#224-jul-5-2022). diff --git a/thirdparty/earcut/include/mapbox/earcut.hpp b/thirdparty/earcut/include/mapbox/earcut.hpp new file mode 100644 index 00000000..b98508f1 --- /dev/null +++ b/thirdparty/earcut/include/mapbox/earcut.hpp @@ -0,0 +1,813 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { + +namespace util { + +template struct nth { + inline static typename std::tuple_element::type + get(const T& t) { return std::get(t); }; +}; + +} + +namespace detail { + +template +class Earcut { +public: + std::vector indices; + std::size_t vertices = 0; + + template + void operator()(const Polygon& points); + +private: + struct Node { + Node(N index, double x_, double y_) : i(index), x(x_), y(y_) {} + Node(const Node&) = delete; + Node& operator=(const Node&) = delete; + Node(Node&&) = delete; + Node& operator=(Node&&) = delete; + + const N i; + const double x; + const double y; + + // previous and next vertice nodes in a polygon ring + Node* prev = nullptr; + Node* next = nullptr; + + // z-order curve value + int32_t z = 0; + + // previous and next nodes in z-order + Node* prevZ = nullptr; + Node* nextZ = nullptr; + + // indicates whether this is a steiner point + bool steiner = false; + }; + + template Node* linkedList(const Ring& points, const bool clockwise); + Node* filterPoints(Node* start, Node* end = nullptr); + void earcutLinked(Node* ear, int pass = 0); + bool isEar(Node* ear); + bool isEarHashed(Node* ear); + Node* cureLocalIntersections(Node* start); + void splitEarcut(Node* start); + template Node* eliminateHoles(const Polygon& points, Node* outerNode); + Node* eliminateHole(Node* hole, Node* outerNode); + Node* findHoleBridge(Node* hole, Node* outerNode); + bool sectorContainsSector(const Node* m, const Node* p); + void indexCurve(Node* start); + Node* sortLinked(Node* list); + int32_t zOrder(const double x_, const double y_); + Node* getLeftmost(Node* start); + bool pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const; + bool isValidDiagonal(Node* a, Node* b); + double area(const Node* p, const Node* q, const Node* r) const; + bool equals(const Node* p1, const Node* p2); + bool intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2); + bool onSegment(const Node* p, const Node* q, const Node* r); + int sign(double val); + bool intersectsPolygon(const Node* a, const Node* b); + bool locallyInside(const Node* a, const Node* b); + bool middleInside(const Node* a, const Node* b); + Node* splitPolygon(Node* a, Node* b); + template Node* insertNode(std::size_t i, const Point& p, Node* last); + void removeNode(Node* p); + + bool hashing; + double minX, maxX; + double minY, maxY; + double inv_size = 0; + + template > + class ObjectPool { + public: + ObjectPool() { } + ObjectPool(std::size_t blockSize_) { + reset(blockSize_); + } + ~ObjectPool() { + clear(); + } + template + T* construct(Args&&... args) { + if (currentIndex >= blockSize) { + currentBlock = alloc_traits::allocate(alloc, blockSize); + allocations.emplace_back(currentBlock); + currentIndex = 0; + } + T* object = ¤tBlock[currentIndex++]; + alloc_traits::construct(alloc, object, std::forward(args)...); + return object; + } + void reset(std::size_t newBlockSize) { + for (auto allocation : allocations) { + alloc_traits::deallocate(alloc, allocation, blockSize); + } + allocations.clear(); + blockSize = std::max(1, newBlockSize); + currentBlock = nullptr; + currentIndex = blockSize; + } + void clear() { reset(blockSize); } + private: + T* currentBlock = nullptr; + std::size_t currentIndex = 1; + std::size_t blockSize = 1; + std::vector allocations; + Alloc alloc; + typedef typename std::allocator_traits alloc_traits; + }; + ObjectPool nodes; +}; + +template template +void Earcut::operator()(const Polygon& points) { + // reset + indices.clear(); + vertices = 0; + + if (points.empty()) return; + + double x; + double y; + int threshold = 80; + std::size_t len = 0; + + for (size_t i = 0; threshold >= 0 && i < points.size(); i++) { + threshold -= static_cast(points[i].size()); + len += points[i].size(); + } + + //estimate size of nodes and indices + nodes.reset(len * 3 / 2); + indices.reserve(len + points[0].size()); + + Node* outerNode = linkedList(points[0], true); + if (!outerNode || outerNode->prev == outerNode->next) return; + + if (points.size() > 1) outerNode = eliminateHoles(points, outerNode); + + // if the shape is not too simple, we'll use z-order curve hash later; calculate polygon bbox + hashing = threshold < 0; + if (hashing) { + Node* p = outerNode->next; + minX = maxX = outerNode->x; + minY = maxY = outerNode->y; + do { + x = p->x; + y = p->y; + minX = std::min(minX, x); + minY = std::min(minY, y); + maxX = std::max(maxX, x); + maxY = std::max(maxY, y); + p = p->next; + } while (p != outerNode); + + // minX, minY and inv_size are later used to transform coords into integers for z-order calculation + inv_size = std::max(maxX - minX, maxY - minY); + inv_size = inv_size != .0 ? (32767. / inv_size) : .0; + } + + earcutLinked(outerNode); + + nodes.clear(); +} + +// create a circular doubly linked list from polygon points in the specified winding order +template template +typename Earcut::Node* +Earcut::linkedList(const Ring& points, const bool clockwise) { + using Point = typename Ring::value_type; + double sum = 0; + const std::size_t len = points.size(); + std::size_t i, j; + Node* last = nullptr; + + // calculate original winding order of a polygon ring + for (i = 0, j = len > 0 ? len - 1 : 0; i < len; j = i++) { + const auto& p1 = points[i]; + const auto& p2 = points[j]; + const double p20 = util::nth<0, Point>::get(p2); + const double p10 = util::nth<0, Point>::get(p1); + const double p11 = util::nth<1, Point>::get(p1); + const double p21 = util::nth<1, Point>::get(p2); + sum += (p20 - p10) * (p11 + p21); + } + + // link points into circular doubly-linked list in the specified winding order + if (clockwise == (sum > 0)) { + for (i = 0; i < len; i++) last = insertNode(vertices + i, points[i], last); + } else { + for (i = len; i-- > 0;) last = insertNode(vertices + i, points[i], last); + } + + if (last && equals(last, last->next)) { + removeNode(last); + last = last->next; + } + + vertices += len; + + return last; +} + +// eliminate colinear or duplicate points +template +typename Earcut::Node* +Earcut::filterPoints(Node* start, Node* end) { + if (!end) end = start; + + Node* p = start; + bool again; + do { + again = false; + + if (!p->steiner && (equals(p, p->next) || area(p->prev, p, p->next) == 0)) { + removeNode(p); + p = end = p->prev; + + if (p == p->next) break; + again = true; + + } else { + p = p->next; + } + } while (again || p != end); + + return end; +} + +// main ear slicing loop which triangulates a polygon (given as a linked list) +template +void Earcut::earcutLinked(Node* ear, int pass) { + if (!ear) return; + + // interlink polygon nodes in z-order + if (!pass && hashing) indexCurve(ear); + + Node* stop = ear; + Node* prev; + Node* next; + + // iterate through ears, slicing them one by one + while (ear->prev != ear->next) { + prev = ear->prev; + next = ear->next; + + if (hashing ? isEarHashed(ear) : isEar(ear)) { + // cut off the triangle + indices.emplace_back(prev->i); + indices.emplace_back(ear->i); + indices.emplace_back(next->i); + + removeNode(ear); + + // skipping the next vertice leads to less sliver triangles + ear = next->next; + stop = next->next; + + continue; + } + + ear = next; + + // if we looped through the whole remaining polygon and can't find any more ears + if (ear == stop) { + // try filtering points and slicing again + if (!pass) earcutLinked(filterPoints(ear), 1); + + // if this didn't work, try curing all small self-intersections locally + else if (pass == 1) { + ear = cureLocalIntersections(filterPoints(ear)); + earcutLinked(ear, 2); + + // as a last resort, try splitting the remaining polygon into two + } else if (pass == 2) splitEarcut(ear); + + break; + } + } +} + +// check whether a polygon node forms a valid ear with adjacent nodes +template +bool Earcut::isEar(Node* ear) { + const Node* a = ear->prev; + const Node* b = ear; + const Node* c = ear->next; + + if (area(a, b, c) >= 0) return false; // reflex, can't be an ear + + // now make sure we don't have other points inside the potential ear + Node* p = ear->next->next; + + while (p != ear->prev) { + if (pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) && + area(p->prev, p, p->next) >= 0) return false; + p = p->next; + } + + return true; +} + +template +bool Earcut::isEarHashed(Node* ear) { + const Node* a = ear->prev; + const Node* b = ear; + const Node* c = ear->next; + + if (area(a, b, c) >= 0) return false; // reflex, can't be an ear + + // triangle bbox; min & max are calculated like this for speed + const double minTX = std::min(a->x, std::min(b->x, c->x)); + const double minTY = std::min(a->y, std::min(b->y, c->y)); + const double maxTX = std::max(a->x, std::max(b->x, c->x)); + const double maxTY = std::max(a->y, std::max(b->y, c->y)); + + // z-order range for the current triangle bbox; + const int32_t minZ = zOrder(minTX, minTY); + const int32_t maxZ = zOrder(maxTX, maxTY); + + // first look for points inside the triangle in increasing z-order + Node* p = ear->nextZ; + + while (p && p->z <= maxZ) { + if (p != ear->prev && p != ear->next && + pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) && + area(p->prev, p, p->next) >= 0) return false; + p = p->nextZ; + } + + // then look for points in decreasing z-order + p = ear->prevZ; + + while (p && p->z >= minZ) { + if (p != ear->prev && p != ear->next && + pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) && + area(p->prev, p, p->next) >= 0) return false; + p = p->prevZ; + } + + return true; +} + +// go through all polygon nodes and cure small local self-intersections +template +typename Earcut::Node* +Earcut::cureLocalIntersections(Node* start) { + Node* p = start; + do { + Node* a = p->prev; + Node* b = p->next->next; + + // a self-intersection where edge (v[i-1],v[i]) intersects (v[i+1],v[i+2]) + if (!equals(a, b) && intersects(a, p, p->next, b) && locallyInside(a, b) && locallyInside(b, a)) { + indices.emplace_back(a->i); + indices.emplace_back(p->i); + indices.emplace_back(b->i); + + // remove two nodes involved + removeNode(p); + removeNode(p->next); + + p = start = b; + } + p = p->next; + } while (p != start); + + return filterPoints(p); +} + +// try splitting polygon into two and triangulate them independently +template +void Earcut::splitEarcut(Node* start) { + // look for a valid diagonal that divides the polygon into two + Node* a = start; + do { + Node* b = a->next->next; + while (b != a->prev) { + if (a->i != b->i && isValidDiagonal(a, b)) { + // split the polygon in two by the diagonal + Node* c = splitPolygon(a, b); + + // filter colinear points around the cuts + a = filterPoints(a, a->next); + c = filterPoints(c, c->next); + + // run earcut on each half + earcutLinked(a); + earcutLinked(c); + return; + } + b = b->next; + } + a = a->next; + } while (a != start); +} + +// link every hole into the outer loop, producing a single-ring polygon without holes +template template +typename Earcut::Node* +Earcut::eliminateHoles(const Polygon& points, Node* outerNode) { + const size_t len = points.size(); + + std::vector queue; + for (size_t i = 1; i < len; i++) { + Node* list = linkedList(points[i], false); + if (list) { + if (list == list->next) list->steiner = true; + queue.push_back(getLeftmost(list)); + } + } + std::sort(queue.begin(), queue.end(), [](const Node* a, const Node* b) { + return a->x < b->x; + }); + + // process holes from left to right + for (size_t i = 0; i < queue.size(); i++) { + outerNode = eliminateHole(queue[i], outerNode); + } + + return outerNode; +} + +// find a bridge between vertices that connects hole with an outer ring and and link it +template +typename Earcut::Node* +Earcut::eliminateHole(Node* hole, Node* outerNode) { + Node* bridge = findHoleBridge(hole, outerNode); + if (!bridge) { + return outerNode; + } + + Node* bridgeReverse = splitPolygon(bridge, hole); + + // filter collinear points around the cuts + filterPoints(bridgeReverse, bridgeReverse->next); + + // Check if input node was removed by the filtering + return filterPoints(bridge, bridge->next); +} + +// David Eberly's algorithm for finding a bridge between hole and outer polygon +template +typename Earcut::Node* +Earcut::findHoleBridge(Node* hole, Node* outerNode) { + Node* p = outerNode; + double hx = hole->x; + double hy = hole->y; + double qx = -std::numeric_limits::infinity(); + Node* m = nullptr; + + // find a segment intersected by a ray from the hole's leftmost Vertex to the left; + // segment's endpoint with lesser x will be potential connection Vertex + do { + if (hy <= p->y && hy >= p->next->y && p->next->y != p->y) { + double x = p->x + (hy - p->y) * (p->next->x - p->x) / (p->next->y - p->y); + if (x <= hx && x > qx) { + qx = x; + m = p->x < p->next->x ? p : p->next; + if (x == hx) return m; // hole touches outer segment; pick leftmost endpoint + } + } + p = p->next; + } while (p != outerNode); + + if (!m) return 0; + + // look for points inside the triangle of hole Vertex, segment intersection and endpoint; + // if there are no points found, we have a valid connection; + // otherwise choose the Vertex of the minimum angle with the ray as connection Vertex + + const Node* stop = m; + double tanMin = std::numeric_limits::infinity(); + double tanCur = 0; + + p = m; + double mx = m->x; + double my = m->y; + + do { + if (hx >= p->x && p->x >= mx && hx != p->x && + pointInTriangle(hy < my ? hx : qx, hy, mx, my, hy < my ? qx : hx, hy, p->x, p->y)) { + + tanCur = std::abs(hy - p->y) / (hx - p->x); // tangential + + if (locallyInside(p, hole) && + (tanCur < tanMin || (tanCur == tanMin && (p->x > m->x || sectorContainsSector(m, p))))) { + m = p; + tanMin = tanCur; + } + } + + p = p->next; + } while (p != stop); + + return m; +} + +// whether sector in vertex m contains sector in vertex p in the same coordinates +template +bool Earcut::sectorContainsSector(const Node* m, const Node* p) { + return area(m->prev, m, p->prev) < 0 && area(p->next, m, m->next) < 0; +} + +// interlink polygon nodes in z-order +template +void Earcut::indexCurve(Node* start) { + assert(start); + Node* p = start; + + do { + p->z = p->z ? p->z : zOrder(p->x, p->y); + p->prevZ = p->prev; + p->nextZ = p->next; + p = p->next; + } while (p != start); + + p->prevZ->nextZ = nullptr; + p->prevZ = nullptr; + + sortLinked(p); +} + +// Simon Tatham's linked list merge sort algorithm +// http://www.chiark.greenend.org.uk/~sgtatham/algorithms/listsort.html +template +typename Earcut::Node* +Earcut::sortLinked(Node* list) { + assert(list); + Node* p; + Node* q; + Node* e; + Node* tail; + int i, numMerges, pSize, qSize; + int inSize = 1; + + for (;;) { + p = list; + list = nullptr; + tail = nullptr; + numMerges = 0; + + while (p) { + numMerges++; + q = p; + pSize = 0; + for (i = 0; i < inSize; i++) { + pSize++; + q = q->nextZ; + if (!q) break; + } + + qSize = inSize; + + while (pSize > 0 || (qSize > 0 && q)) { + + if (pSize == 0) { + e = q; + q = q->nextZ; + qSize--; + } else if (qSize == 0 || !q) { + e = p; + p = p->nextZ; + pSize--; + } else if (p->z <= q->z) { + e = p; + p = p->nextZ; + pSize--; + } else { + e = q; + q = q->nextZ; + qSize--; + } + + if (tail) tail->nextZ = e; + else list = e; + + e->prevZ = tail; + tail = e; + } + + p = q; + } + + tail->nextZ = nullptr; + + if (numMerges <= 1) return list; + + inSize *= 2; + } +} + +// z-order of a Vertex given coords and size of the data bounding box +template +int32_t Earcut::zOrder(const double x_, const double y_) { + // coords are transformed into non-negative 15-bit integer range + int32_t x = static_cast((x_ - minX) * inv_size); + int32_t y = static_cast((y_ - minY) * inv_size); + + x = (x | (x << 8)) & 0x00FF00FF; + x = (x | (x << 4)) & 0x0F0F0F0F; + x = (x | (x << 2)) & 0x33333333; + x = (x | (x << 1)) & 0x55555555; + + y = (y | (y << 8)) & 0x00FF00FF; + y = (y | (y << 4)) & 0x0F0F0F0F; + y = (y | (y << 2)) & 0x33333333; + y = (y | (y << 1)) & 0x55555555; + + return x | (y << 1); +} + +// find the leftmost node of a polygon ring +template +typename Earcut::Node* +Earcut::getLeftmost(Node* start) { + Node* p = start; + Node* leftmost = start; + do { + if (p->x < leftmost->x || (p->x == leftmost->x && p->y < leftmost->y)) + leftmost = p; + p = p->next; + } while (p != start); + + return leftmost; +} + +// check if a point lies within a convex triangle +template +bool Earcut::pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const { + return (cx - px) * (ay - py) >= (ax - px) * (cy - py) && + (ax - px) * (by - py) >= (bx - px) * (ay - py) && + (bx - px) * (cy - py) >= (cx - px) * (by - py); +} + +// check if a diagonal between two polygon nodes is valid (lies in polygon interior) +template +bool Earcut::isValidDiagonal(Node* a, Node* b) { + return a->next->i != b->i && a->prev->i != b->i && !intersectsPolygon(a, b) && // dones't intersect other edges + ((locallyInside(a, b) && locallyInside(b, a) && middleInside(a, b) && // locally visible + (area(a->prev, a, b->prev) != 0.0 || area(a, b->prev, b) != 0.0)) || // does not create opposite-facing sectors + (equals(a, b) && area(a->prev, a, a->next) > 0 && area(b->prev, b, b->next) > 0)); // special zero-length case +} + +// signed area of a triangle +template +double Earcut::area(const Node* p, const Node* q, const Node* r) const { + return (q->y - p->y) * (r->x - q->x) - (q->x - p->x) * (r->y - q->y); +} + +// check if two points are equal +template +bool Earcut::equals(const Node* p1, const Node* p2) { + return p1->x == p2->x && p1->y == p2->y; +} + +// check if two segments intersect +template +bool Earcut::intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2) { + int o1 = sign(area(p1, q1, p2)); + int o2 = sign(area(p1, q1, q2)); + int o3 = sign(area(p2, q2, p1)); + int o4 = sign(area(p2, q2, q1)); + + if (o1 != o2 && o3 != o4) return true; // general case + + if (o1 == 0 && onSegment(p1, p2, q1)) return true; // p1, q1 and p2 are collinear and p2 lies on p1q1 + if (o2 == 0 && onSegment(p1, q2, q1)) return true; // p1, q1 and q2 are collinear and q2 lies on p1q1 + if (o3 == 0 && onSegment(p2, p1, q2)) return true; // p2, q2 and p1 are collinear and p1 lies on p2q2 + if (o4 == 0 && onSegment(p2, q1, q2)) return true; // p2, q2 and q1 are collinear and q1 lies on p2q2 + + return false; +} + +// for collinear points p, q, r, check if point q lies on segment pr +template +bool Earcut::onSegment(const Node* p, const Node* q, const Node* r) { + return q->x <= std::max(p->x, r->x) && + q->x >= std::min(p->x, r->x) && + q->y <= std::max(p->y, r->y) && + q->y >= std::min(p->y, r->y); +} + +template +int Earcut::sign(double val) { + return (0.0 < val) - (val < 0.0); +} + +// check if a polygon diagonal intersects any polygon segments +template +bool Earcut::intersectsPolygon(const Node* a, const Node* b) { + const Node* p = a; + do { + if (p->i != a->i && p->next->i != a->i && p->i != b->i && p->next->i != b->i && + intersects(p, p->next, a, b)) return true; + p = p->next; + } while (p != a); + + return false; +} + +// check if a polygon diagonal is locally inside the polygon +template +bool Earcut::locallyInside(const Node* a, const Node* b) { + return area(a->prev, a, a->next) < 0 ? + area(a, b, a->next) >= 0 && area(a, a->prev, b) >= 0 : + area(a, b, a->prev) < 0 || area(a, a->next, b) < 0; +} + +// check if the middle Vertex of a polygon diagonal is inside the polygon +template +bool Earcut::middleInside(const Node* a, const Node* b) { + const Node* p = a; + bool inside = false; + double px = (a->x + b->x) / 2; + double py = (a->y + b->y) / 2; + do { + if (((p->y > py) != (p->next->y > py)) && p->next->y != p->y && + (px < (p->next->x - p->x) * (py - p->y) / (p->next->y - p->y) + p->x)) + inside = !inside; + p = p->next; + } while (p != a); + + return inside; +} + +// link two polygon vertices with a bridge; if the vertices belong to the same ring, it splits +// polygon into two; if one belongs to the outer ring and another to a hole, it merges it into a +// single ring +template +typename Earcut::Node* +Earcut::splitPolygon(Node* a, Node* b) { + Node* a2 = nodes.construct(a->i, a->x, a->y); + Node* b2 = nodes.construct(b->i, b->x, b->y); + Node* an = a->next; + Node* bp = b->prev; + + a->next = b; + b->prev = a; + + a2->next = an; + an->prev = a2; + + b2->next = a2; + a2->prev = b2; + + bp->next = b2; + b2->prev = bp; + + return b2; +} + +// create a node and util::optionally link it with previous one (in a circular doubly linked list) +template template +typename Earcut::Node* +Earcut::insertNode(std::size_t i, const Point& pt, Node* last) { + Node* p = nodes.construct(static_cast(i), util::nth<0, Point>::get(pt), util::nth<1, Point>::get(pt)); + + if (!last) { + p->prev = p; + p->next = p; + + } else { + assert(last); + p->next = last->next; + p->prev = last; + last->next->prev = p; + last->next = p; + } + return p; +} + +template +void Earcut::removeNode(Node* p) { + p->next->prev = p->prev; + p->prev->next = p->next; + + if (p->prevZ) p->prevZ->nextZ = p->nextZ; + if (p->nextZ) p->nextZ->prevZ = p->prevZ; +} +} + +template +std::vector earcut(const Polygon& poly) { + mapbox::detail::Earcut earcut; + earcut(poly); + return std::move(earcut.indices); +} +}