polygon.cpp | Explore in Territory

// Copyright 2021 The Manifold Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//      http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "manifold/polygon.h"

#include <functional>
#include <map>
#include <set>

#include "./collider.h"
#include "./parallel.h"
#include "./utils.h"
#include "manifold/optional_assert.h"

namespace {
usingnamespacemanifold;

static ExecutionParams params;

constexpr double kBest = …;

// it seems that MSVC cannot optimize la::determinant(mat2(a, b))
constexpr double determinant2x2(vec2 a, vec2 b) { … }

#ifdef MANIFOLD_DEBUG
struct PolyEdge {
  int startVert, endVert;
};

std::vector<PolyEdge> Polygons2Edges(const PolygonsIdx &polys) {
  std::vector<PolyEdge> halfedges;
  for (const auto &poly : polys) {
    for (size_t i = 1; i < poly.size(); ++i) {
      halfedges.push_back({poly[i - 1].idx, poly[i].idx});
    }
    halfedges.push_back({poly.back().idx, poly[0].idx});
  }
  return halfedges;
}

std::vector<PolyEdge> Triangles2Edges(const std::vector<ivec3> &triangles) {
  std::vector<PolyEdge> halfedges;
  halfedges.reserve(triangles.size() * 3);
  for (const ivec3 &tri : triangles) {
    halfedges.push_back({tri[0], tri[1]});
    halfedges.push_back({tri[1], tri[2]});
    halfedges.push_back({tri[2], tri[0]});
  }
  return halfedges;
}

void CheckTopology(const std::vector<PolyEdge> &halfedges) {
  DEBUG_ASSERT(halfedges.size() % 2 == 0, topologyErr,
               "Odd number of halfedges.");
  size_t n_edges = halfedges.size() / 2;
  std::vector<PolyEdge> forward(halfedges.size()), backward(halfedges.size());

  auto end = std::copy_if(halfedges.begin(), halfedges.end(), forward.begin(),
                          [](PolyEdge e) { return e.endVert > e.startVert; });
  DEBUG_ASSERT(
      static_cast<size_t>(std::distance(forward.begin(), end)) == n_edges,
      topologyErr, "Half of halfedges should be forward.");
  forward.resize(n_edges);

  end = std::copy_if(halfedges.begin(), halfedges.end(), backward.begin(),
                     [](PolyEdge e) { return e.endVert < e.startVert; });
  DEBUG_ASSERT(
      static_cast<size_t>(std::distance(backward.begin(), end)) == n_edges,
      topologyErr, "Half of halfedges should be backward.");
  backward.resize(n_edges);

  std::for_each(backward.begin(), backward.end(),
                [](PolyEdge &e) { std::swap(e.startVert, e.endVert); });
  auto cmp = [](const PolyEdge &a, const PolyEdge &b) {
    return a.startVert < b.startVert ||
           (a.startVert == b.startVert && a.endVert < b.endVert);
  };
  std::stable_sort(forward.begin(), forward.end(), cmp);
  std::stable_sort(backward.begin(), backward.end(), cmp);
  for (size_t i = 0; i < n_edges; ++i) {
    DEBUG_ASSERT(forward[i].startVert == backward[i].startVert &&
                     forward[i].endVert == backward[i].endVert,
                 topologyErr, "Not manifold.");
  }
}

void CheckTopology(const std::vector<ivec3> &triangles,
                   const PolygonsIdx &polys) {
  std::vector<PolyEdge> halfedges = Triangles2Edges(triangles);
  std::vector<PolyEdge> openEdges = Polygons2Edges(polys);
  for (PolyEdge e : openEdges) {
    halfedges.push_back({e.endVert, e.startVert});
  }
  CheckTopology(halfedges);
}

void CheckGeometry(const std::vector<ivec3> &triangles,
                   const PolygonsIdx &polys, double epsilon) {
  std::unordered_map<int, vec2> vertPos;
  for (const auto &poly : polys) {
    for (size_t i = 0; i < poly.size(); ++i) {
      vertPos[poly[i].idx] = poly[i].pos;
    }
  }
  DEBUG_ASSERT(std::all_of(triangles.begin(), triangles.end(),
                           [&vertPos, epsilon](const ivec3 &tri) {
                             return CCW(vertPos[tri[0]], vertPos[tri[1]],
                                        vertPos[tri[2]], epsilon) >= 0;
                           }),
               geometryErr, "triangulation is not entirely CCW!");
}

void Dump(const PolygonsIdx &polys, double epsilon) {
  std::cout << "Polygon 0 " << epsilon << " " << polys.size() << std::endl;
  for (auto poly : polys) {
    std::cout << poly.size() << std::endl;
    for (auto v : poly) {
      std::cout << v.pos.x << " " << v.pos.y << std::endl;
    }
  }
  std::cout << "# ... " << std::endl;
  for (auto poly : polys) {
    std::cout << "show(array([" << std::endl;
    for (auto v : poly) {
      std::cout << "  [" << v.pos.x << ", " << v.pos.y << "]," << std::endl;
    }
    std::cout << "]))" << std::endl;
  }
}

void PrintFailure(const std::exception &e, const PolygonsIdx &polys,
                  std::vector<ivec3> &triangles, double epsilon) {
  std::cout << "-----------------------------------" << std::endl;
  std::cout << "Triangulation failed! Precision = " << epsilon << std::endl;
  std::cout << e.what() << std::endl;
  if (triangles.size() > 1000 && !PolygonParams().verbose) {
    std::cout << "Output truncated due to producing " << triangles.size()
              << " triangles." << std::endl;
    return;
  }
  Dump(polys, epsilon);
  std::cout << "produced this triangulation:" << std::endl;
  for (size_t j = 0; j < triangles.size(); ++j) {
    std::cout << triangles[j][0] << ", " << triangles[j][1] << ", "
              << triangles[j][2] << std::endl;
  }
}

#define PRINT …
#else
#define PRINT(msg) …
#endif

/**
 * Tests if the input polygons are convex by searching for any reflex vertices.
 * Exactly colinear edges and zero-length edges are treated conservatively as
 * reflex. Does not check for overlaps.
 */
bool IsConvex(const PolygonsIdx &polys, double epsilon) { … }

/**
 * Triangulates a set of convex polygons by alternating instead of a fan, to
 * avoid creating high-degree vertices.
 */
std::vector<ivec3> TriangulateConvex(const PolygonsIdx &polys) { … }

/**
 * Ear-clipping triangulator based on David Eberly's approach from Geometric
 * Tools, but adjusted to handle epsilon-valid polygons, and including a
 * fallback that ensures a manifold triangulation even for overlapping polygons.
 * This is an O(n^2) algorithm, but hopefully this is not a big problem as the
 * number of edges in a given polygon is generally much less than the number of
 * triangles in a mesh, and relatively few faces even need triangulation.
 *
 * The main adjustments for robustness involve clipping the sharpest ears first
 * (a known technique to get higher triangle quality), and doing an exhaustive
 * search to determine ear convexity exactly if the first geometric result is
 * within epsilon.
 */

class EarClip { … };
}  // namespace

namespace manifold {

/**
 * @brief Triangulates a set of &epsilon;-valid polygons. If the input is not
 * &epsilon;-valid, the triangulation may overlap, but will always return a
 * manifold result that matches the input edge directions.
 *
 * @param polys The set of polygons, wound CCW and representing multiple
 * polygons and/or holes. These have 2D-projected positions as well as
 * references back to the original vertices.
 * @param epsilon The value of &epsilon;, bounding the uncertainty of the
 * input.
 * @return std::vector<ivec3> The triangles, referencing the original
 * vertex indicies.
 */
std::vector<ivec3> TriangulateIdx(const PolygonsIdx &polys, double epsilon) { … }

/**
 * @brief Triangulates a set of &epsilon;-valid polygons. If the input is not
 * &epsilon;-valid, the triangulation may overlap, but will always return a
 * manifold result that matches the input edge directions.
 *
 * @param polygons The set of polygons, wound CCW and representing multiple
 * polygons and/or holes.
 * @param epsilon The value of &epsilon;, bounding the uncertainty of the
 * input.
 * @return std::vector<ivec3> The triangles, referencing the original
 * polygon points in order.
 */
std::vector<ivec3> Triangulate(const Polygons &polygons, double epsilon) { … }

ExecutionParams &PolygonParams() { … }

}  // namespace manifold
godot/thirdparty/manifold/src/polygon.cpp
