// Copyright 2023 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/cross_section.h" #include "../utils.h" #include "clipper2/clipper.core.h" #include "clipper2/clipper.h" #include "clipper2/clipper.offset.h" namespace C2 = Clipper2Lib; using namespace manifold; namespace manifold { struct PathImpl { PathImpl(const C2::PathsD paths_) : paths_(paths_) {} operator const C2::PathsD&() const { return paths_; } const C2::PathsD paths_; }; } // namespace manifold namespace { const int precision_ = 8; C2::ClipType cliptype_of_op(OpType op) { C2::ClipType ct = C2::ClipType::Union; switch (op) { case OpType::Add: break; case OpType::Subtract: ct = C2::ClipType::Difference; break; case OpType::Intersect: ct = C2::ClipType::Intersection; break; }; return ct; } C2::FillRule fr(CrossSection::FillRule fillrule) { C2::FillRule fr = C2::FillRule::EvenOdd; switch (fillrule) { case CrossSection::FillRule::EvenOdd: break; case CrossSection::FillRule::NonZero: fr = C2::FillRule::NonZero; break; case CrossSection::FillRule::Positive: fr = C2::FillRule::Positive; break; case CrossSection::FillRule::Negative: fr = C2::FillRule::Negative; break; }; return fr; } C2::JoinType jt(CrossSection::JoinType jointype) { C2::JoinType jt = C2::JoinType::Square; switch (jointype) { case CrossSection::JoinType::Square: break; case CrossSection::JoinType::Round: jt = C2::JoinType::Round; break; case CrossSection::JoinType::Miter: jt = C2::JoinType::Miter; break; case CrossSection::JoinType::Bevel: jt = C2::JoinType::Bevel; break; }; return jt; } vec2 v2_of_pd(const C2::PointD p) { return {p.x, p.y}; } C2::PointD v2_to_pd(const vec2 v) { return C2::PointD(v.x, v.y); } C2::PathD pathd_of_contour(const SimplePolygon& ctr) { auto p = C2::PathD(); p.reserve(ctr.size()); for (auto v : ctr) { p.push_back(v2_to_pd(v)); } return p; } C2::PathsD transform(const C2::PathsD ps, const mat2x3 m) { const bool invert = la::determinant(mat2(m)) < 0; auto transformed = C2::PathsD(); transformed.reserve(ps.size()); for (auto path : ps) { auto sz = path.size(); auto s = C2::PathD(sz); for (size_t i = 0; i < sz; ++i) { auto idx = invert ? sz - 1 - i : i; s[idx] = v2_to_pd(m * vec3(path[i].x, path[i].y, 1)); } transformed.push_back(s); } return transformed; } std::shared_ptr shared_paths(const C2::PathsD& ps) { return std::make_shared(ps); } // forward declaration for mutual recursion void decompose_hole(const C2::PolyTreeD* outline, std::vector& polys, C2::PathsD& poly, size_t n_holes, size_t j); void decompose_outline(const C2::PolyTreeD* tree, std::vector& polys, size_t i) { auto n_outlines = tree->Count(); if (i < n_outlines) { auto outline = tree->Child(i); auto n_holes = outline->Count(); auto poly = C2::PathsD(n_holes + 1); poly[0] = outline->Polygon(); decompose_hole(outline, polys, poly, n_holes, 0); polys.push_back(poly); if (i < n_outlines - 1) { decompose_outline(tree, polys, i + 1); } } } void decompose_hole(const C2::PolyTreeD* outline, std::vector& polys, C2::PathsD& poly, size_t n_holes, size_t j) { if (j < n_holes) { auto child = outline->Child(j); decompose_outline(child, polys, 0); poly[j + 1] = child->Polygon(); decompose_hole(outline, polys, poly, n_holes, j + 1); } } void flatten(const C2::PolyTreeD* tree, C2::PathsD& polys, size_t i) { auto n_outlines = tree->Count(); if (i < n_outlines) { auto outline = tree->Child(i); flatten(outline, polys, 0); polys.push_back(outline->Polygon()); if (i < n_outlines - 1) { flatten(tree, polys, i + 1); } } } bool V2Lesser(vec2 a, vec2 b) { if (a.x == b.x) return a.y < b.y; return a.x < b.x; } void HullBacktrack(const vec2& pt, std::vector& stack) { auto sz = stack.size(); while (sz >= 2 && CCW(stack[sz - 2], stack[sz - 1], pt, 0.0) <= 0.0) { stack.pop_back(); sz = stack.size(); } } // Based on method described here: // https://www.hackerearth.com/practice/math/geometry/line-sweep-technique/tutorial/ // Changed to follow: // https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain // This is the same algorithm (Andrew, also called Montone Chain). C2::PathD HullImpl(SimplePolygon& pts) { size_t len = pts.size(); if (len < 3) return C2::PathD(); // not enough points to create a polygon std::sort(pts.begin(), pts.end(), V2Lesser); auto lower = std::vector{}; for (auto& pt : pts) { HullBacktrack(pt, lower); lower.push_back(pt); } auto upper = std::vector{}; for (auto pt_iter = pts.rbegin(); pt_iter != pts.rend(); pt_iter++) { HullBacktrack(*pt_iter, upper); upper.push_back(*pt_iter); } upper.pop_back(); lower.pop_back(); auto path = C2::PathD(); path.reserve(lower.size() + upper.size()); for (const auto& l : lower) path.push_back(v2_to_pd(l)); for (const auto& u : upper) path.push_back(v2_to_pd(u)); return path; } } // namespace namespace manifold { /** * The default constructor is an empty cross-section (containing no contours). */ CrossSection::CrossSection() { paths_ = std::make_shared(C2::PathsD()); } CrossSection::~CrossSection() = default; CrossSection::CrossSection(CrossSection&&) noexcept = default; CrossSection& CrossSection::operator=(CrossSection&&) noexcept = default; /** * The copy constructor avoids copying the underlying paths vector (sharing * with its parent via shared_ptr), however subsequent transformations, and * their application will not be shared. It is generally recommended to avoid * this, opting instead to simply create CrossSections with the available * const methods. */ CrossSection::CrossSection(const CrossSection& other) { paths_ = other.paths_; transform_ = other.transform_; } CrossSection& CrossSection::operator=(const CrossSection& other) { if (this != &other) { paths_ = other.paths_; transform_ = other.transform_; } return *this; }; // Private, skips unioning. CrossSection::CrossSection(std::shared_ptr ps) { paths_ = ps; } /** * Create a 2d cross-section from a single contour. A boolean union operation * (with Positive filling rule by default) is performed to ensure the * resulting CrossSection is free of self-intersections. * * @param contour A closed path outlining the desired cross-section. * @param fillrule The filling rule used to interpret polygon sub-regions * created by self-intersections in contour. */ CrossSection::CrossSection(const SimplePolygon& contour, FillRule fillrule) { auto ps = C2::PathsD{(pathd_of_contour(contour))}; paths_ = shared_paths(C2::Union(ps, fr(fillrule), precision_)); } /** * Create a 2d cross-section from a set of contours (complex polygons). A * boolean union operation (with Positive filling rule by default) is * performed to combine overlapping polygons and ensure the resulting * CrossSection is free of intersections. * * @param contours A set of closed paths describing zero or more complex * polygons. * @param fillrule The filling rule used to interpret polygon sub-regions in * contours. */ CrossSection::CrossSection(const Polygons& contours, FillRule fillrule) { auto ps = C2::PathsD(); ps.reserve(contours.size()); for (auto ctr : contours) { ps.push_back(pathd_of_contour(ctr)); } paths_ = shared_paths(C2::Union(ps, fr(fillrule), precision_)); } /** * Create a 2d cross-section from an axis-aligned rectangle (bounding box). * * @param rect An axis-aligned rectangular bounding box. */ CrossSection::CrossSection(const Rect& rect) { C2::PathD p(4); p[0] = C2::PointD(rect.min.x, rect.min.y); p[1] = C2::PointD(rect.max.x, rect.min.y); p[2] = C2::PointD(rect.max.x, rect.max.y); p[3] = C2::PointD(rect.min.x, rect.max.y); paths_ = shared_paths(C2::PathsD{p}); } // Private // All access to paths_ should be done through the GetPaths() method, which // applies the accumulated transform_ std::shared_ptr CrossSection::GetPaths() const { if (transform_ == mat2x3(la::identity)) { return paths_; } paths_ = shared_paths(::transform(paths_->paths_, transform_)); transform_ = mat2x3(la::identity); return paths_; } /** * Constructs a square with the given XY dimensions. By default it is * positioned in the first quadrant, touching the origin. If any dimensions in * size are negative, or if all are zero, an empty Manifold will be returned. * * @param size The X, and Y dimensions of the square. * @param center Set to true to shift the center to the origin. */ CrossSection CrossSection::Square(const vec2 size, bool center) { if (size.x < 0.0 || size.y < 0.0 || la::length(size) == 0.0) { return CrossSection(); } auto p = C2::PathD(4); if (center) { const auto w = size.x / 2; const auto h = size.y / 2; p[0] = C2::PointD(w, h); p[1] = C2::PointD(-w, h); p[2] = C2::PointD(-w, -h); p[3] = C2::PointD(w, -h); } else { const double x = size.x; const double y = size.y; p[0] = C2::PointD(0.0, 0.0); p[1] = C2::PointD(x, 0.0); p[2] = C2::PointD(x, y); p[3] = C2::PointD(0.0, y); } return CrossSection(shared_paths(C2::PathsD{p})); } /** * Constructs a circle of a given radius. * * @param radius Radius of the circle. Must be positive. * @param circularSegments Number of segments along its diameter. Default is * calculated by the static Quality defaults according to the radius. */ CrossSection CrossSection::Circle(double radius, int circularSegments) { if (radius <= 0.0) { return CrossSection(); } int n = circularSegments > 2 ? circularSegments : Quality::GetCircularSegments(radius); double dPhi = 360.0 / n; auto circle = C2::PathD(n); for (int i = 0; i < n; ++i) { circle[i] = C2::PointD(radius * cosd(dPhi * i), radius * sind(dPhi * i)); } return CrossSection(shared_paths(C2::PathsD{circle})); } /** * Perform the given boolean operation between this and another CrossSection. */ CrossSection CrossSection::Boolean(const CrossSection& second, OpType op) const { auto ct = cliptype_of_op(op); auto res = C2::BooleanOp(ct, C2::FillRule::Positive, GetPaths()->paths_, second.GetPaths()->paths_, precision_); return CrossSection(shared_paths(res)); } /** * Perform the given boolean operation on a list of CrossSections. In case of * Subtract, all CrossSections in the tail are differenced from the head. */ CrossSection CrossSection::BatchBoolean( const std::vector& crossSections, OpType op) { if (crossSections.size() == 0) return CrossSection(); else if (crossSections.size() == 1) return crossSections[0]; auto subjs = crossSections[0].GetPaths(); if (op == OpType::Intersect) { auto res = subjs->paths_; for (size_t i = 1; i < crossSections.size(); ++i) { res = C2::BooleanOp(C2::ClipType::Intersection, C2::FillRule::Positive, res, crossSections[i].GetPaths()->paths_, precision_); } return CrossSection(shared_paths(res)); } int n_clips = 0; for (size_t i = 1; i < crossSections.size(); ++i) { n_clips += crossSections[i].GetPaths()->paths_.size(); } auto clips = C2::PathsD(); clips.reserve(n_clips); for (size_t i = 1; i < crossSections.size(); ++i) { auto ps = crossSections[i].GetPaths(); clips.insert(clips.end(), ps->paths_.begin(), ps->paths_.end()); } auto ct = cliptype_of_op(op); auto res = C2::BooleanOp(ct, C2::FillRule::Positive, subjs->paths_, clips, precision_); return CrossSection(shared_paths(res)); } /** * Compute the boolean union between two cross-sections. */ CrossSection CrossSection::operator+(const CrossSection& Q) const { return Boolean(Q, OpType::Add); } /** * Compute the boolean union between two cross-sections, assigning the result * to the first. */ CrossSection& CrossSection::operator+=(const CrossSection& Q) { *this = *this + Q; return *this; } /** * Compute the boolean difference of a (clip) cross-section from another * (subject). */ CrossSection CrossSection::operator-(const CrossSection& Q) const { return Boolean(Q, OpType::Subtract); } /** * Compute the boolean difference of a (clip) cross-section from a another * (subject), assigning the result to the subject. */ CrossSection& CrossSection::operator-=(const CrossSection& Q) { *this = *this - Q; return *this; } /** * Compute the boolean intersection between two cross-sections. */ CrossSection CrossSection::operator^(const CrossSection& Q) const { return Boolean(Q, OpType::Intersect); } /** * Compute the boolean intersection between two cross-sections, assigning the * result to the first. */ CrossSection& CrossSection::operator^=(const CrossSection& Q) { *this = *this ^ Q; return *this; } /** * Construct a CrossSection from a vector of other CrossSections (batch * boolean union). */ CrossSection CrossSection::Compose( const std::vector& crossSections) { return BatchBoolean(crossSections, OpType::Add); } /** * This operation returns a vector of CrossSections that are topologically * disconnected, each containing one outline contour with zero or more * holes. */ std::vector CrossSection::Decompose() const { if (NumContour() < 2) { return std::vector{CrossSection(*this)}; } C2::PolyTreeD tree; C2::BooleanOp(C2::ClipType::Union, C2::FillRule::Positive, GetPaths()->paths_, C2::PathsD(), tree, precision_); auto polys = std::vector(); decompose_outline(&tree, polys, 0); auto comps = std::vector(); comps.reserve(polys.size()); // reverse the stack while wrapping for (auto poly = polys.rbegin(); poly != polys.rend(); ++poly) comps.emplace_back(CrossSection(shared_paths(*poly))); return comps; } /** * Move this CrossSection in space. This operation can be chained. Transforms * are combined and applied lazily. * * @param v The vector to add to every vertex. */ CrossSection CrossSection::Translate(const vec2 v) const { mat2x3 m({1.0, 0.0}, // {0.0, 1.0}, // {v.x, v.y}); return Transform(m); } /** * Applies a (Z-axis) rotation to the CrossSection, in degrees. This operation * can be chained. Transforms are combined and applied lazily. * * @param degrees degrees about the Z-axis to rotate. */ CrossSection CrossSection::Rotate(double degrees) const { auto s = sind(degrees); auto c = cosd(degrees); mat2x3 m({c, s}, // {-s, c}, // {0.0, 0.0}); return Transform(m); } /** * Scale this CrossSection in space. This operation can be chained. Transforms * are combined and applied lazily. * * @param scale The vector to multiply every vertex by per component. */ CrossSection CrossSection::Scale(const vec2 scale) const { mat2x3 m({scale.x, 0.0}, // {0.0, scale.y}, // {0.0, 0.0}); return Transform(m); } /** * Mirror this CrossSection over the arbitrary axis whose normal is described by * the unit form of the given vector. If the length of the vector is zero, an * empty CrossSection is returned. This operation can be chained. Transforms are * combined and applied lazily. * * @param ax the axis to be mirrored over */ CrossSection CrossSection::Mirror(const vec2 ax) const { if (la::length(ax) == 0.) { return CrossSection(); } auto n = la::normalize(ax); auto m = mat2x3(mat2(la::identity) - 2.0 * la::outerprod(n, n), vec2(0.0)); return Transform(m); } /** * Transform this CrossSection in space. The first two columns form a 2x2 * matrix transform and the last is a translation vector. This operation can * be chained. Transforms are combined and applied lazily. * * @param m The affine transform matrix to apply to all the vertices. */ CrossSection CrossSection::Transform(const mat2x3& m) const { auto transformed = CrossSection(); transformed.transform_ = m * Mat3(transform_); transformed.paths_ = paths_; return transformed; } /** * Move the vertices of this CrossSection (creating a new one) according to * any arbitrary input function, followed by a union operation (with a * Positive fill rule) that ensures any introduced intersections are not * included in the result. * * @param warpFunc A function that modifies a given vertex position. */ CrossSection CrossSection::Warp(std::function warpFunc) const { return WarpBatch([&warpFunc](VecView vecs) { for (vec2& p : vecs) { warpFunc(p); } }); } /** * Same as CrossSection::Warp but calls warpFunc with * a VecView which is roughly equivalent to std::span * pointing to all vec2 elements to be modified in-place * * @param warpFunc A function that modifies multiple vertex positions. */ CrossSection CrossSection::WarpBatch( std::function)> warpFunc) const { std::vector tmp_verts; C2::PathsD paths = GetPaths()->paths_; // deep copy for (C2::PathD const& path : paths) { for (C2::PointD const& p : path) { tmp_verts.push_back(v2_of_pd(p)); } } warpFunc(VecView(tmp_verts.data(), tmp_verts.size())); auto cursor = tmp_verts.begin(); for (C2::PathD& path : paths) { for (C2::PointD& p : path) { p = v2_to_pd(*cursor); ++cursor; } } return CrossSection( shared_paths(C2::Union(paths, C2::FillRule::Positive, precision_))); } /** * Remove vertices from the contours in this CrossSection that are less than * the specified distance epsilon from an imaginary line that passes through * its two adjacent vertices. Near duplicate vertices and collinear points * will be removed at lower epsilons, with elimination of line segments * becoming increasingly aggressive with larger epsilons. * * It is recommended to apply this function following Offset, in order to * clean up any spurious tiny line segments introduced that do not improve * quality in any meaningful way. This is particularly important if further * offseting operations are to be performed, which would compound the issue. */ CrossSection CrossSection::Simplify(double epsilon) const { C2::PolyTreeD tree; C2::BooleanOp(C2::ClipType::Union, C2::FillRule::Positive, GetPaths()->paths_, C2::PathsD(), tree, precision_); C2::PathsD polys; flatten(&tree, polys, 0); // Filter out contours less than epsilon wide. C2::PathsD filtered; for (C2::PathD poly : polys) { auto area = C2::Area(poly); Rect box; for (auto vert : poly) { box.Union(vec2(vert.x, vert.y)); } vec2 size = box.Size(); if (std::abs(area) > std::max(size.x, size.y) * epsilon) { filtered.push_back(poly); } } auto ps = SimplifyPaths(filtered, epsilon, true); return CrossSection(shared_paths(ps)); } /** * Inflate the contours in CrossSection by the specified delta, handling * corners according to the given JoinType. * * @param delta Positive deltas will cause the expansion of outlining contours * to expand, and retraction of inner (hole) contours. Negative deltas will * have the opposite effect. * @param jointype The join type specifying the treatment of contour joins * (corners). Defaults to Round. * @param miter_limit The maximum distance in multiples of delta that vertices * can be offset from their original positions with before squaring is * applied, when the join type is Miter (default is 2, which is the * minimum allowed). See the [Clipper2 * MiterLimit](http://www.angusj.com/clipper2/Docs/Units/Clipper.Offset/Classes/ClipperOffset/Properties/MiterLimit.htm) * page for a visual example. * @param circularSegments Number of segments per 360 degrees of * JoinType::Round corners (roughly, the number of vertices that * will be added to each contour). Default is calculated by the static Quality * defaults according to the radius. */ CrossSection CrossSection::Offset(double delta, JoinType jointype, double miter_limit, int circularSegments) const { double arc_tol = 0.; if (jointype == JoinType::Round) { int n = circularSegments > 2 ? circularSegments : Quality::GetCircularSegments(delta); // This calculates tolerance as a function of circular segments and delta // (radius) in order to get back the same number of segments in Clipper2: // steps_per_360 = PI / acos(1 - arc_tol / abs_delta) const double abs_delta = std::fabs(delta); arc_tol = (std::cos(Clipper2Lib::PI / n) - 1) * -abs_delta; } auto ps = C2::InflatePaths(GetPaths()->paths_, delta, jt(jointype), C2::EndType::Polygon, miter_limit, precision_, arc_tol); return CrossSection(shared_paths(ps)); } /** * Compute the convex hull enveloping a set of cross-sections. * * @param crossSections A vector of cross-sections over which to compute a * convex hull. */ CrossSection CrossSection::Hull( const std::vector& crossSections) { int n = 0; for (auto cs : crossSections) n += cs.NumVert(); SimplePolygon pts; pts.reserve(n); for (auto cs : crossSections) { auto paths = cs.GetPaths()->paths_; for (auto path : paths) { for (auto p : path) { pts.push_back(v2_of_pd(p)); } } } return CrossSection(shared_paths(C2::PathsD{HullImpl(pts)})); } /** * Compute the convex hull of this cross-section. */ CrossSection CrossSection::Hull() const { return Hull(std::vector{*this}); } /** * Compute the convex hull of a set of points. If the given points are fewer * than 3, an empty CrossSection will be returned. * * @param pts A vector of 2-dimensional points over which to compute a convex * hull. */ CrossSection CrossSection::Hull(SimplePolygon pts) { return CrossSection(shared_paths(C2::PathsD{HullImpl(pts)})); } /** * Compute the convex hull of a set of points/polygons. If the given points are * fewer than 3, an empty CrossSection will be returned. * * @param polys A vector of vectors of 2-dimensional points over which to * compute a convex hull. */ CrossSection CrossSection::Hull(const Polygons polys) { SimplePolygon pts; for (auto poly : polys) { for (auto p : poly) { pts.push_back(p); } } return Hull(pts); } /** * Return the total area covered by complex polygons making up the * CrossSection. */ double CrossSection::Area() const { return C2::Area(GetPaths()->paths_); } /** * Return the number of vertices in the CrossSection. */ size_t CrossSection::NumVert() const { size_t n = 0; auto paths = GetPaths()->paths_; for (auto p : paths) { n += p.size(); } return n; } /** * Return the number of contours (both outer and inner paths) in the * CrossSection. */ size_t CrossSection::NumContour() const { return GetPaths()->paths_.size(); } /** * Does the CrossSection contain any contours? */ bool CrossSection::IsEmpty() const { return GetPaths()->paths_.empty(); } /** * Returns the axis-aligned bounding rectangle of all the CrossSections' * vertices. */ Rect CrossSection::Bounds() const { auto r = C2::GetBounds(GetPaths()->paths_); return Rect({r.left, r.bottom}, {r.right, r.top}); } /** * Return the contours of this CrossSection as a Polygons. */ Polygons CrossSection::ToPolygons() const { auto polys = Polygons(); auto paths = GetPaths()->paths_; polys.reserve(paths.size()); for (auto p : paths) { auto sp = SimplePolygon(); sp.reserve(p.size()); for (auto v : p) { sp.push_back({v.x, v.y}); } polys.push_back(sp); } return polys; } } // namespace manifold