File indexing completed on 2024-05-12 03:52:08

 /*******************************************************************************
 * Author    :  Angus Johnson                                                   *
 * Date      :  24 September 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::DoBevel(Group& group, const Path64& path, size_t j, size_t k)
 {
     PointD pt1, pt2;
     if (j == k)
     {
         double abs_delta = std::abs(group_delta_);
         pt1 = PointD(path[j].x - abs_delta * norms[j].x, path[j].y - abs_delta * norms[j].y);
         pt2 = PointD(path[j].x + abs_delta * norms[j].x, path[j].y + abs_delta * norms[j].y);
     } 
     else
     {
         pt1 = PointD(path[j].x + group_delta_ * norms[k].x, path[j].y + group_delta_ * norms[k].y);
         pt2 = PointD(path[j].x + group_delta_ * norms[j].x, path[j].y + group_delta_ * norms[j].y);
     }
     group.path.push_back(Point64(pt1));
     group.path.push_back(Point64(pt2));
 }
 
 void ClipperOffset::DoSquare(Group& group, const Path64& path, size_t j, size_t k)
 {
     PointD vec;
     if (j == k) 
         vec = PointD(norms[j].y, -norms[j].x);
     else
         vec = GetAvgUnitVector(
             PointD(-norms[k].y, norms[k].x),
             PointD(norms[j].y, -norms[j].x));
 
     double abs_delta = std::abs(group_delta_);
 
     // now offset the original vertex delta units along unit vector
     PointD ptQ = PointD(path[j]);
     ptQ = TranslatePoint(ptQ, abs_delta * vec.x, abs_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)
 {
     if (deltaCallback64_) {
         // when deltaCallback64_ is assigned, group_delta_ won't be constant, 
         // so we'll need to do the following calculations for *every* vertex.
         double abs_delta = std::fabs(group_delta_);
         double arcTol = (arc_tolerance_ > floating_point_tolerance ?
             std::min(abs_delta, arc_tolerance_) :
             std::log10(2 + abs_delta) * default_arc_tolerance);
         double steps_per_360 = std::min(PI / std::acos(1 - arcTol / abs_delta), abs_delta * PI);
         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);
     }
 
     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
     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 (deltaCallback64_) {
         group_delta_ = deltaCallback64_(path, norms, j, k);
         if (group.is_reversed) group_delta_ = -group_delta_;
     }
     if (std::fabs(group_delta_) <= floating_point_tolerance)
     {
         group.path.push_back(path[j]);
         return;
     }
 
     if (cos_a > -0.99 && (sin_a * group_delta_ < 0)) // test for concavity first (#593)
     {
         // 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 (cos_a > 0.999) // almost straight - less than 2.5 degree (#424, #526) 
     {
         DoMiter(group, path, j, k, 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);
     }
     else if (cos_a > 0.99 || join_type_ == JoinType::Bevel) 
         // ie > 2.5 deg (see above) but less than ~8 deg ( acos(0.99) )
         DoBevel(group, path, j, k);
     else if (join_type_ == JoinType::Round)
         DoRound(group, path, j, k, std::atan2(sin_a, cos_a));
     else
         DoSquare(group, path, j, k);
 }
 
 void ClipperOffset::OffsetPolygon(Group& group, Path64& path)
 {
     // when the path is contracting, make sure 
     // there is sufficient space to do so.                //#593
     // nb: this will have a small impact on performance
     double a = Area(path);
     // contracting when orientation is opposite offset direction
     if ((a < 0) != (group_delta_ < 0)) 
     {
         Rect64 rec = GetBounds(path);
         double offsetMinDim = std::fabs(group_delta_) * 2;
         if (offsetMinDim > rec.Width() || offsetMinDim > rec.Height()) return;
     }
 
     for (Path64::size_type j = 0, k = path.size() -1; j < path.size(); k = j, ++j)
         OffsetPoint(group, path, j, k);
     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
     if (deltaCallback64_) group_delta_ = deltaCallback64_(path, norms, 0, 0);
     
     if (std::fabs(group_delta_) <= floating_point_tolerance)
         group.path.push_back(path[0]);
     else
     {
         switch (end_type_)
         {
         case EndType::Butt:
             DoBevel(group, path, 0, 0);
             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 j = 1, k = 0; j < highI; k = j, ++j)
         OffsetPoint(group, path, j, 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
     if (deltaCallback64_)
         group_delta_ = deltaCallback64_(path, norms, highI, highI);
 
     if (std::fabs(group_delta_) <= floating_point_tolerance)
         group.path.push_back(path[highI]);
     else
     {
         switch (end_type_)
         {
         case EndType::Butt:
             DoBevel(group, path, highI, highI);
             break;
         case EndType::Round:
             DoRound(group, path, highI, highI, PI);
             break;
         default:
             DoSquare(group, path, highI, highI);
             break;
         }
     }
 
     for (size_t j = highI, k = 0; j > 0; k = j, --j)
         OffsetPoint(group, path, j, 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;
     }
 
     double abs_delta = std::fabs(group_delta_);
     // do range checking
     if (!IsSafeOffset(r, abs_delta))
     {
         DoError(range_error_i);
         error_code_ |= range_error_i;
         return;
     }
 
     join_type_  = group.join_type;
     end_type_ = group.end_type;
 
     if (!deltaCallback64_ && 
         (group.join_type == JoinType::Round || group.end_type == EndType::Round))
     {
         //calculate a sensible number of steps (for 360 deg for the given offset)
         // 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_delta, arc_tolerance_) :
             std::log10(2 + abs_delta) * default_arc_tolerance);
 
         double steps_per_360 = std::min(PI / std::acos(1 - arcTol / abs_delta), abs_delta * PI);
         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::iterator path_iter;
     for(path_iter = group.paths_in.begin(); path_iter != group.paths_in.end(); ++path_iter)
     {
         Path64 &path = *path_iter;
         StripDuplicates(path, 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_delta;
                 int steps = static_cast<int>(std::ceil(steps_per_rad_ * 2 * PI)); //#617
                 group.path = Ellipse(path[0], radius, radius, steps);
 #ifdef USINGZ
                 for (auto& p : group.path) p.z = path[0].z;
 #endif
             }
             else
             {
                 int d = (int)std::ceil(abs_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);
 }
 
 void ClipperOffset::Execute(DeltaCallback64 delta_cb, Paths64& paths)
 {
     deltaCallback64_ = delta_cb;
     Execute(1.0, paths);
 }
 
 } // namespace