diff options
Diffstat (limited to 'thirdparty/clipper2/src/clipper.offset.cpp')
| -rw-r--r-- | thirdparty/clipper2/src/clipper.offset.cpp | 618 |
1 files changed, 618 insertions, 0 deletions
diff --git a/thirdparty/clipper2/src/clipper.offset.cpp b/thirdparty/clipper2/src/clipper.offset.cpp new file mode 100644 index 0000000000..78cd82376a --- /dev/null +++ b/thirdparty/clipper2/src/clipper.offset.cpp @@ -0,0 +1,618 @@ +/******************************************************************************* +* Author : Angus Johnson * +* Date : 22 March 2023 * +* Website : http://www.angusj.com * +* Copyright : Angus Johnson 2010-2023 * +* Purpose : Path Offset (Inflate/Shrink) * +* License : http://www.boost.org/LICENSE_1_0.txt * +*******************************************************************************/ + +#include <cmath> +#include "clipper2/clipper.h" +#include "clipper2/clipper.offset.h" + +namespace Clipper2Lib { + +const double default_arc_tolerance = 0.25; +const double floating_point_tolerance = 1e-12; + +//------------------------------------------------------------------------------ +// Miscellaneous methods +//------------------------------------------------------------------------------ + +void GetBoundsAndLowestPolyIdx(const Paths64& paths, Rect64& r, int & idx) +{ + idx = -1; + r = MaxInvalidRect64; + int64_t lpx = 0; + for (int i = 0; i < static_cast<int>(paths.size()); ++i) + for (const Point64& p : paths[i]) + { + if (p.y >= r.bottom) + { + if (p.y > r.bottom || p.x < lpx) + { + idx = i; + lpx = p.x; + r.bottom = p.y; + } + } + else if (p.y < r.top) r.top = p.y; + if (p.x > r.right) r.right = p.x; + else if (p.x < r.left) r.left = p.x; + } + //if (idx < 0) r = Rect64(0, 0, 0, 0); + //if (r.top == INT64_MIN) r.bottom = r.top; + //if (r.left == INT64_MIN) r.left = r.right; +} + +bool IsSafeOffset(const Rect64& r, double abs_delta) +{ + return r.left > min_coord + abs_delta && + r.right < max_coord - abs_delta && + r.top > min_coord + abs_delta && + r.bottom < max_coord - abs_delta; +} + +PointD GetUnitNormal(const Point64& pt1, const Point64& pt2) +{ + double dx, dy, inverse_hypot; + if (pt1 == pt2) return PointD(0.0, 0.0); + dx = static_cast<double>(pt2.x - pt1.x); + dy = static_cast<double>(pt2.y - pt1.y); + inverse_hypot = 1.0 / hypot(dx, dy); + dx *= inverse_hypot; + dy *= inverse_hypot; + return PointD(dy, -dx); +} + +inline bool AlmostZero(double value, double epsilon = 0.001) +{ + return std::fabs(value) < epsilon; +} + +inline double Hypot(double x, double y) +{ + //see https://stackoverflow.com/a/32436148/359538 + return std::sqrt(x * x + y * y); +} + +inline PointD NormalizeVector(const PointD& vec) +{ + + double h = Hypot(vec.x, vec.y); + if (AlmostZero(h)) return PointD(0,0); + double inverseHypot = 1 / h; + return PointD(vec.x * inverseHypot, vec.y * inverseHypot); +} + +inline PointD GetAvgUnitVector(const PointD& vec1, const PointD& vec2) +{ + return NormalizeVector(PointD(vec1.x + vec2.x, vec1.y + vec2.y)); +} + +inline bool IsClosedPath(EndType et) +{ + return et == EndType::Polygon || et == EndType::Joined; +} + +inline Point64 GetPerpendic(const Point64& pt, const PointD& norm, double delta) +{ +#ifdef USINGZ + return Point64(pt.x + norm.x * delta, pt.y + norm.y * delta, pt.z); +#else + return Point64(pt.x + norm.x * delta, pt.y + norm.y * delta); +#endif +} + +inline PointD GetPerpendicD(const Point64& pt, const PointD& norm, double delta) +{ +#ifdef USINGZ + return PointD(pt.x + norm.x * delta, pt.y + norm.y * delta, pt.z); +#else + return PointD(pt.x + norm.x * delta, pt.y + norm.y * delta); +#endif +} + +inline void NegatePath(PathD& path) +{ + for (PointD& pt : path) + { + pt.x = -pt.x; + pt.y = -pt.y; +#ifdef USINGZ + pt.z = pt.z; +#endif + } +} + +//------------------------------------------------------------------------------ +// ClipperOffset methods +//------------------------------------------------------------------------------ + +void ClipperOffset::AddPath(const Path64& path, JoinType jt_, EndType et_) +{ + Paths64 paths; + paths.push_back(path); + AddPaths(paths, jt_, et_); +} + +void ClipperOffset::AddPaths(const Paths64 &paths, JoinType jt_, EndType et_) +{ + if (paths.size() == 0) return; + groups_.push_back(Group(paths, jt_, et_)); +} + +void ClipperOffset::BuildNormals(const Path64& path) +{ + norms.clear(); + norms.reserve(path.size()); + if (path.size() == 0) return; + Path64::const_iterator path_iter, path_last_iter = --path.cend(); + for (path_iter = path.cbegin(); path_iter != path_last_iter; ++path_iter) + norms.push_back(GetUnitNormal(*path_iter,*(path_iter +1))); + norms.push_back(GetUnitNormal(*path_last_iter, *(path.cbegin()))); +} + +inline PointD TranslatePoint(const PointD& pt, double dx, double dy) +{ +#ifdef USINGZ + return PointD(pt.x + dx, pt.y + dy, pt.z); +#else + return PointD(pt.x + dx, pt.y + dy); +#endif +} + +inline PointD ReflectPoint(const PointD& pt, const PointD& pivot) +{ +#ifdef USINGZ + return PointD(pivot.x + (pivot.x - pt.x), pivot.y + (pivot.y - pt.y), pt.z); +#else + return PointD(pivot.x + (pivot.x - pt.x), pivot.y + (pivot.y - pt.y)); +#endif +} + +PointD IntersectPoint(const PointD& pt1a, const PointD& pt1b, + const PointD& pt2a, const PointD& pt2b) +{ + if (pt1a.x == pt1b.x) //vertical + { + if (pt2a.x == pt2b.x) return PointD(0, 0); + + double m2 = (pt2b.y - pt2a.y) / (pt2b.x - pt2a.x); + double b2 = pt2a.y - m2 * pt2a.x; + return PointD(pt1a.x, m2 * pt1a.x + b2); + } + else if (pt2a.x == pt2b.x) //vertical + { + double m1 = (pt1b.y - pt1a.y) / (pt1b.x - pt1a.x); + double b1 = pt1a.y - m1 * pt1a.x; + return PointD(pt2a.x, m1 * pt2a.x + b1); + } + else + { + double m1 = (pt1b.y - pt1a.y) / (pt1b.x - pt1a.x); + double b1 = pt1a.y - m1 * pt1a.x; + double m2 = (pt2b.y - pt2a.y) / (pt2b.x - pt2a.x); + double b2 = pt2a.y - m2 * pt2a.x; + if (m1 == m2) return PointD(0, 0); + double x = (b2 - b1) / (m1 - m2); + return PointD(x, m1 * x + b1); + } +} + +void ClipperOffset::DoSquare(Group& group, const Path64& path, size_t j, size_t k) +{ + PointD vec; + if (j == k) + vec = PointD(norms[0].y, -norms[0].x); + else + vec = GetAvgUnitVector( + PointD(-norms[k].y, norms[k].x), + PointD(norms[j].y, -norms[j].x)); + + // now offset the original vertex delta units along unit vector + PointD ptQ = PointD(path[j]); + ptQ = TranslatePoint(ptQ, abs_group_delta_ * vec.x, abs_group_delta_ * vec.y); + // get perpendicular vertices + PointD pt1 = TranslatePoint(ptQ, group_delta_ * vec.y, group_delta_ * -vec.x); + PointD pt2 = TranslatePoint(ptQ, group_delta_ * -vec.y, group_delta_ * vec.x); + // get 2 vertices along one edge offset + PointD pt3 = GetPerpendicD(path[k], norms[k], group_delta_); + if (j == k) + { + PointD pt4 = PointD(pt3.x + vec.x * group_delta_, pt3.y + vec.y * group_delta_); + PointD pt = IntersectPoint(pt1, pt2, pt3, pt4); +#ifdef USINGZ + pt.z = ptQ.z; +#endif + //get the second intersect point through reflecion + group.path.push_back(Point64(ReflectPoint(pt, ptQ))); + group.path.push_back(Point64(pt)); + } + else + { + PointD pt4 = GetPerpendicD(path[j], norms[k], group_delta_); + PointD pt = IntersectPoint(pt1, pt2, pt3, pt4); +#ifdef USINGZ + pt.z = ptQ.z; +#endif + group.path.push_back(Point64(pt)); + //get the second intersect point through reflecion + group.path.push_back(Point64(ReflectPoint(pt, ptQ))); + } +} + +void ClipperOffset::DoMiter(Group& group, const Path64& path, size_t j, size_t k, double cos_a) +{ + double q = group_delta_ / (cos_a + 1); +#ifdef USINGZ + group.path.push_back(Point64( + path[j].x + (norms[k].x + norms[j].x) * q, + path[j].y + (norms[k].y + norms[j].y) * q, + path[j].z)); +#else + group.path.push_back(Point64( + path[j].x + (norms[k].x + norms[j].x) * q, + path[j].y + (norms[k].y + norms[j].y) * q)); +#endif +} + +void ClipperOffset::DoRound(Group& group, const Path64& path, size_t j, size_t k, double angle) +{ + Point64 pt = path[j]; + PointD offsetVec = PointD(norms[k].x * group_delta_, norms[k].y * group_delta_); + + if (j == k) offsetVec.Negate(); +#ifdef USINGZ + group.path.push_back(Point64(pt.x + offsetVec.x, pt.y + offsetVec.y, pt.z)); +#else + group.path.push_back(Point64(pt.x + offsetVec.x, pt.y + offsetVec.y)); +#endif + if (angle > -PI + 0.01) // avoid 180deg concave + { + int steps = static_cast<int>(std::ceil(steps_per_rad_ * std::abs(angle))); // #448, #456 + for (int i = 1; i < steps; ++i) // ie 1 less than steps + { + offsetVec = PointD(offsetVec.x * step_cos_ - step_sin_ * offsetVec.y, + offsetVec.x * step_sin_ + offsetVec.y * step_cos_); +#ifdef USINGZ + group.path.push_back(Point64(pt.x + offsetVec.x, pt.y + offsetVec.y, pt.z)); +#else + group.path.push_back(Point64(pt.x + offsetVec.x, pt.y + offsetVec.y)); +#endif + + } + } + group.path.push_back(GetPerpendic(path[j], norms[j], group_delta_)); +} + +void ClipperOffset::OffsetPoint(Group& group, Path64& path, size_t j, size_t& k) +{ + // Let A = change in angle where edges join + // A == 0: ie no change in angle (flat join) + // A == PI: edges 'spike' + // sin(A) < 0: right turning + // cos(A) < 0: change in angle is more than 90 degree + + if (path[j] == path[k]) { k = j; return; } + + double sin_a = CrossProduct(norms[j], norms[k]); + double cos_a = DotProduct(norms[j], norms[k]); + if (sin_a > 1.0) sin_a = 1.0; + else if (sin_a < -1.0) sin_a = -1.0; + + if (cos_a > 0.99) // almost straight - less than 8 degrees + { + group.path.push_back(GetPerpendic(path[j], norms[k], group_delta_)); + if (cos_a < 0.9998) // greater than 1 degree (#424) + group.path.push_back(GetPerpendic(path[j], norms[j], group_delta_)); // (#418) + } + else if (cos_a > -0.99 && (sin_a * group_delta_ < 0)) + { + // is concave + group.path.push_back(GetPerpendic(path[j], norms[k], group_delta_)); + // this extra point is the only (simple) way to ensure that + // path reversals are fully cleaned with the trailing clipper + group.path.push_back(path[j]); // (#405) + group.path.push_back(GetPerpendic(path[j], norms[j], group_delta_)); + } + else if (join_type_ == JoinType::Round) + DoRound(group, path, j, k, std::atan2(sin_a, cos_a)); + else if (join_type_ == JoinType::Miter) + { + // miter unless the angle is so acute the miter would exceeds ML + if (cos_a > temp_lim_ - 1) DoMiter(group, path, j, k, cos_a); + else DoSquare(group, path, j, k); + } + // don't bother squaring angles that deviate < ~20 degrees because + // squaring will be indistinguishable from mitering and just be a lot slower + else if (cos_a > 0.9) + DoMiter(group, path, j, k, cos_a); + else + DoSquare(group, path, j, k); + + k = j; +} + +void ClipperOffset::OffsetPolygon(Group& group, Path64& path) +{ + for (Path64::size_type i = 0, j = path.size() -1; i < path.size(); j = i, ++i) + OffsetPoint(group, path, i, j); + group.paths_out.push_back(group.path); +} + +void ClipperOffset::OffsetOpenJoined(Group& group, Path64& path) +{ + OffsetPolygon(group, path); + std::reverse(path.begin(), path.end()); + + //rebuild normals // BuildNormals(path); + std::reverse(norms.begin(), norms.end()); + norms.push_back(norms[0]); + norms.erase(norms.begin()); + NegatePath(norms); + + group.path.clear(); + OffsetPolygon(group, path); +} + +void ClipperOffset::OffsetOpenPath(Group& group, Path64& path) +{ + // do the line start cap + switch (end_type_) + { + case EndType::Butt: +#ifdef USINGZ + group.path.push_back(Point64( + path[0].x - norms[0].x * group_delta_, + path[0].y - norms[0].y * group_delta_, + path[0].z)); +#else + group.path.push_back(Point64( + path[0].x - norms[0].x * group_delta_, + path[0].y - norms[0].y * group_delta_)); +#endif + group.path.push_back(GetPerpendic(path[0], norms[0], group_delta_)); + break; + case EndType::Round: + DoRound(group, path, 0,0, PI); + break; + default: + DoSquare(group, path, 0, 0); + break; + } + + size_t highI = path.size() - 1; + + // offset the left side going forward + for (Path64::size_type i = 1, k = 0; i < highI; ++i) + OffsetPoint(group, path, i, k); + + // reverse normals + for (size_t i = highI; i > 0; --i) + norms[i] = PointD(-norms[i - 1].x, -norms[i - 1].y); + norms[0] = norms[highI]; + + // do the line end cap + switch (end_type_) + { + case EndType::Butt: +#ifdef USINGZ + group.path.push_back(Point64( + path[highI].x - norms[highI].x * group_delta_, + path[highI].y - norms[highI].y * group_delta_, + path[highI].z)); +#else + group.path.push_back(Point64( + path[highI].x - norms[highI].x * group_delta_, + path[highI].y - norms[highI].y * group_delta_)); +#endif + group.path.push_back(GetPerpendic(path[highI], norms[highI], group_delta_)); + break; + case EndType::Round: + DoRound(group, path, highI, highI, PI); + break; + default: + DoSquare(group, path, highI, highI); + break; + } + + for (size_t i = highI, k = 0; i > 0; --i) + OffsetPoint(group, path, i, k); + group.paths_out.push_back(group.path); +} + +void ClipperOffset::DoGroupOffset(Group& group) +{ + Rect64 r; + int idx = -1; + //the lowermost polygon must be an outer polygon. So we can use that as the + //designated orientation for outer polygons (needed for tidy-up clipping) + GetBoundsAndLowestPolyIdx(group.paths_in, r, idx); + if (idx < 0) return; + + if (group.end_type == EndType::Polygon) + { + double area = Area(group.paths_in[idx]); + //if (area == 0) return; // probably unhelpful (#430) + group.is_reversed = (area < 0); + if (group.is_reversed) group_delta_ = -delta_; + else group_delta_ = delta_; + } + else + { + group.is_reversed = false; + group_delta_ = std::abs(delta_) * 0.5; + } + abs_group_delta_ = std::fabs(group_delta_); + + // do range checking + if (!IsSafeOffset(r, abs_group_delta_)) + { + DoError(range_error_i); + error_code_ |= range_error_i; + return; + } + + join_type_ = group.join_type; + end_type_ = group.end_type; + + //calculate a sensible number of steps (for 360 deg for the given offset + if (group.join_type == JoinType::Round || group.end_type == EndType::Round) + { + // arcTol - when arc_tolerance_ is undefined (0), the amount of + // curve imprecision that's allowed is based on the size of the + // offset (delta). Obviously very large offsets will almost always + // require much less precision. See also offset_triginometry2.svg + double arcTol = (arc_tolerance_ > floating_point_tolerance ? + std::min(abs_group_delta_, arc_tolerance_) : + std::log10(2 + abs_group_delta_) * default_arc_tolerance); + double steps_per_360 = PI / std::acos(1 - arcTol / abs_group_delta_); + if (steps_per_360 > abs_group_delta_ * PI) + steps_per_360 = abs_group_delta_ * PI; //ie avoids excessive precision + + step_sin_ = std::sin(2 * PI / steps_per_360); + step_cos_ = std::cos(2 * PI / steps_per_360); + if (group_delta_ < 0.0) step_sin_ = -step_sin_; + steps_per_rad_ = steps_per_360 / (2 *PI); + } + + bool is_joined = + (end_type_ == EndType::Polygon) || + (end_type_ == EndType::Joined); + Paths64::const_iterator path_iter; + for(path_iter = group.paths_in.cbegin(); path_iter != group.paths_in.cend(); ++path_iter) + { + Path64 path = StripDuplicates(*path_iter, is_joined); + Path64::size_type cnt = path.size(); + if (cnt == 0 || ((cnt < 3) && group.end_type == EndType::Polygon)) + continue; + + group.path.clear(); + if (cnt == 1) // single point - only valid with open paths + { + if (group_delta_ < 1) continue; + //single vertex so build a circle or square ... + if (group.join_type == JoinType::Round) + { + double radius = abs_group_delta_; + group.path = Ellipse(path[0], radius, radius); +#ifdef USINGZ + for (auto& p : group.path) p.z = path[0].z; +#endif + } + else + { + int d = (int)std::ceil(abs_group_delta_); + r = Rect64(path[0].x - d, path[0].y - d, path[0].x + d, path[0].y + d); + group.path = r.AsPath(); +#ifdef USINGZ + for (auto& p : group.path) p.z = path[0].z; +#endif + } + group.paths_out.push_back(group.path); + } + else + { + if ((cnt == 2) && (group.end_type == EndType::Joined)) + { + if (group.join_type == JoinType::Round) + end_type_ = EndType::Round; + else + end_type_ = EndType::Square; + } + + BuildNormals(path); + if (end_type_ == EndType::Polygon) OffsetPolygon(group, path); + else if (end_type_ == EndType::Joined) OffsetOpenJoined(group, path); + else OffsetOpenPath(group, path); + } + } + solution.reserve(solution.size() + group.paths_out.size()); + copy(group.paths_out.begin(), group.paths_out.end(), back_inserter(solution)); + group.paths_out.clear(); +} + +void ClipperOffset::ExecuteInternal(double delta) +{ + error_code_ = 0; + solution.clear(); + if (groups_.size() == 0) return; + + if (std::abs(delta) < 0.5) + { + for (const Group& group : groups_) + { + solution.reserve(solution.size() + group.paths_in.size()); + copy(group.paths_in.begin(), group.paths_in.end(), back_inserter(solution)); + } + } + else + { + temp_lim_ = (miter_limit_ <= 1) ? + 2.0 : + 2.0 / (miter_limit_ * miter_limit_); + + delta_ = delta; + std::vector<Group>::iterator git; + for (git = groups_.begin(); git != groups_.end(); ++git) + { + DoGroupOffset(*git); + if (!error_code_) continue; // all OK + solution.clear(); + } + } +} + +void ClipperOffset::Execute(double delta, Paths64& paths) +{ + paths.clear(); + + ExecuteInternal(delta); + if (!solution.size()) return; + + paths = solution; + //clean up self-intersections ... + Clipper64 c; + c.PreserveCollinear = false; + //the solution should retain the orientation of the input + c.ReverseSolution = reverse_solution_ != groups_[0].is_reversed; +#ifdef USINGZ + if (zCallback64_) { + c.SetZCallback(zCallback64_); + } +#endif + c.AddSubject(solution); + if (groups_[0].is_reversed) + c.Execute(ClipType::Union, FillRule::Negative, paths); + else + c.Execute(ClipType::Union, FillRule::Positive, paths); +} + + +void ClipperOffset::Execute(double delta, PolyTree64& polytree) +{ + polytree.Clear(); + + ExecuteInternal(delta); + if (!solution.size()) return; + + //clean up self-intersections ... + Clipper64 c; + c.PreserveCollinear = false; + //the solution should retain the orientation of the input + c.ReverseSolution = reverse_solution_ != groups_[0].is_reversed; +#ifdef USINGZ + if (zCallback64_) { + c.SetZCallback(zCallback64_); + } +#endif + c.AddSubject(solution); + if (groups_[0].is_reversed) + c.Execute(ClipType::Union, FillRule::Negative, polytree); + else + c.Execute(ClipType::Union, FillRule::Positive, polytree); +} + +} // namespace |