| File: | librecad/src/lib/engine/document/container/lc_loopclassifier.cpp |
| Warning: | line 558, column 10 Value stored to 'yAtPhi' during its initialization is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | /* |
| 2 | ********************************************************************************** |
| 3 | ** |
| 4 | ** This file was created for the LibreCAD project (librecad.org), a 2D CAD program. |
| 5 | ** |
| 6 | ** Copyright (C) 2026 librecad (www.librecad.org) |
| 7 | ** |
| 8 | ** This program is free software; you can redistribute it and/or |
| 9 | ** modify it under the terms of the GNU General Public License |
| 10 | ** as published by the Free Software Foundation; either version 2 |
| 11 | ** of the License, or (at your option) any later version. |
| 12 | ** |
| 13 | ** This program is distributed in the hope that it will be useful, |
| 14 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 15 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 16 | ** GNU General Public License for more details. |
| 17 | ** |
| 18 | ********************************************************************************** |
| 19 | */ |
| 20 | |
| 21 | #include "lc_loopclassifier.h" |
| 22 | |
| 23 | #include <algorithm> |
| 24 | #include <array> |
| 25 | #include <cfloat> |
| 26 | #include <cmath> |
| 27 | #include <cstdint> |
| 28 | #include <limits> |
| 29 | #include <memory> |
| 30 | #include <numeric> |
| 31 | #include <unordered_map> |
| 32 | #include <vector> |
| 33 | |
| 34 | #include "lc_containertraverser.h" |
| 35 | #include "lc_hyperbola.h" |
| 36 | #include "lc_parabola.h" |
| 37 | #include "lc_rect.h" |
| 38 | #include "lc_splinepoints.h" |
| 39 | #include "rs.h" |
| 40 | #include "rs_arc.h" |
| 41 | #include "rs_circle.h" |
| 42 | #include "rs_debug.h" |
| 43 | #include "rs_ellipse.h" |
| 44 | #include "rs_entity.h" |
| 45 | #include "rs_entitycontainer.h" |
| 46 | #include "rs_line.h" |
| 47 | #include "rs_math.h" |
| 48 | #include "rs_point.h" |
| 49 | #include "rs_vector.h" |
| 50 | |
| 51 | namespace LC_LoopUtils { |
| 52 | |
| 53 | namespace { |
| 54 | |
| 55 | constexpr double kU = DBL_EPSILON2.2204460492503131e-16 / 2.0; // 2^-53 |
| 56 | constexpr double kPi = M_PI3.14159265358979323846; |
| 57 | constexpr double kTwoPi = 2.0 * M_PI3.14159265358979323846; |
| 58 | constexpr double kHalfPi = M_PI3.14159265358979323846 / 2.0; |
| 59 | |
| 60 | // An extremum nearer than this to an end of an arc gets no split. A split |
| 61 | // nearer than the round-off of a join's angle (a few 1e-15) would only add a |
| 62 | // piece whose ends differ by rounding. Beyond sqrt(u) = 1.05e-8 one is |
| 63 | // needed: up to there sin() of an arc's angle is exactly 1 at an extremum, so |
| 64 | // the end vertex's y is the extremum's y to the bit, but further out the |
| 65 | // vertex can sit a rounding step short of it, and without a split no piece |
| 66 | // covers the y values in between. The margin is well above the first bound |
| 67 | // and 11 times below the second. Elliptic arcs use the same value. |
| 68 | constexpr double kSplitMargin = 1.0 / 1073741824.0; // 2^-30 |
| 69 | |
| 70 | // An invalid vector holds (0, 0), so its coordinates alone would pass as a point. |
| 71 | bool isFiniteV(const RS_Vector& v) { |
| 72 | return v.valid && std::isfinite(v.x) && std::isfinite(v.y); |
| 73 | } |
| 74 | |
| 75 | /** Normalise angle to [0, 2*pi). */ |
| 76 | double normalise2Pi(double a) { |
| 77 | double r = std::fmod(a, kTwoPi); |
| 78 | if (r < 0) r += kTwoPi; |
| 79 | return r; |
| 80 | } |
| 81 | |
| 82 | /** sign(x) with sign(0) = 1 (used in stable quadratic formula). */ |
| 83 | double signPositive(double x) { |
| 84 | return x < 0.0 ? -1.0 : 1.0; |
| 85 | } |
| 86 | |
| 87 | /** Pre-filter: is the bbox of entity e within epsOn of p? */ |
| 88 | bool bboxContainsPoint(const RS_Entity& e, const RS_Vector& p, double epsOn) { |
| 89 | const RS_Vector mn = e.getMin(); |
| 90 | const RS_Vector mx = e.getMax(); |
| 91 | return p.x >= mn.x - epsOn && p.x <= mx.x + epsOn && |
| 92 | p.y >= mn.y - epsOn && p.y <= mx.y + epsOn; |
| 93 | } |
| 94 | |
| 95 | } // namespace |
| 96 | |
| 97 | double hatchTolerance(double diagonal, double maxAbsCoordinate) { |
| 98 | const double term1 = 64.0 * kU * maxAbsCoordinate; |
| 99 | const double term2 = 1.0e-9 * diagonal; |
| 100 | const double floor = std::min(1.0e-8, 1.0e-6 * diagonal); |
| 101 | return std::max({term1, term2, floor}); |
| 102 | } |
| 103 | |
| 104 | // -------------------------------------------------------------------------- |
| 105 | // MonoPiece::xAt |
| 106 | // -------------------------------------------------------------------------- |
| 107 | |
| 108 | double MonoPiece::xAt(double y) const { |
| 109 | const double yc = std::clamp(y, ylo, yhi); |
| 110 | if (yc == ylo) return a.x; |
| 111 | if (yc == yhi) return b.x; |
| 112 | |
| 113 | switch (kind) { |
| 114 | case Kind::Line: { |
| 115 | const double dy = yhi - ylo; |
| 116 | if (dy == 0.0) return a.x; |
| 117 | const double t = (yc - ylo) / dy; |
| 118 | return a.x + t * (b.x - a.x); |
| 119 | } |
| 120 | case Kind::Arc: { |
| 121 | const double dy = yc - cy; |
| 122 | const double rr = r * r - dy * dy; |
| 123 | const double root = rr > 0.0 ? std::sqrt(rr) : 0.0; |
| 124 | // xAt = cx + s * (r - dy^2 / (r + sqrt(r^2 - dy^2))) stable form |
| 125 | if (r + root == 0.0) return cx + s * r; |
| 126 | return cx + s * (r - dy * dy / (r + root)); |
| 127 | } |
| 128 | case Kind::Ellipse: { |
| 129 | if (R == 0.0) return a.x; |
| 130 | double v = (yc - cy) / R; |
| 131 | v = std::clamp(v, -1.0, 1.0); |
| 132 | // half == 0: t = alpha + 2*pi*m + acos(v); half == 1: t = alpha + 2*pi*m + 2*pi - acos(v) |
| 133 | // We know the piece is monotone and tLo/tHi bracket t for ylo/yhi. |
| 134 | // We can derive which multiple m to use from tLo. |
| 135 | const double base = (half == 0) ? std::acos(v) : (kTwoPi - std::acos(v)); |
| 136 | // tLo has the form alpha + 2*pi*m + base_at_ylo. We pick m such that |
| 137 | // t is within [min(tLo, tHi), max(tLo, tHi)]. |
| 138 | double t = majorAngle + base; |
| 139 | const double tMin = std::min(tLo, tHi); |
| 140 | const double tMax = std::max(tLo, tHi); |
| 141 | // adjust by 2*pi*m |
| 142 | while (t < tMin - 1e-9) t += kTwoPi; |
| 143 | while (t > tMax + 1e-9) t -= kTwoPi; |
| 144 | t = std::clamp(t, tMin, tMax); |
| 145 | return cx + mx * std::cos(t) + nx * std::sin(t); |
| 146 | } |
| 147 | case Kind::Quadratic: { |
| 148 | // qy(w) = qa*w^2 + qb*w + qc0 = y with qa = y0 - 2y1 + y2, qb = 2(y1-y0), qc0 = y0 |
| 149 | // Rewrite as qa*w^2 + qb*w + (qc0 - y) = 0, i.e. coefficients qa, qb, (qc0 - y) |
| 150 | const double qa = qy0 - 2.0 * qy1 + qy2; |
| 151 | const double qb = 2.0 * (qy1 - qy0); |
| 152 | const double qc = qy0 - yc; |
| 153 | double w; |
| 154 | if (std::abs(qa) < 1e-14 * (std::abs(qb) + std::abs(qc) + 1.0)) { |
| 155 | // linear |
| 156 | if (qb == 0.0) w = wLo; |
| 157 | else w = -qc / qb; |
| 158 | } else { |
| 159 | const double disc = std::max(0.0, qb * qb - 4.0 * qa * qc); |
| 160 | const double q = -(qb + signPositive(qb) * std::sqrt(disc)) / 2.0; |
| 161 | // roots q/qa and qc/q; pick the one inside [wLo, wHi] |
| 162 | const double r1 = q / qa; |
| 163 | const double r2 = (q != 0.0) ? (qc / q) : r1; |
| 164 | const double wmin = std::min(wLo, wHi); |
| 165 | const double wmax = std::max(wLo, wHi); |
| 166 | const bool r1in = (r1 >= wmin - 1e-9 && r1 <= wmax + 1e-9); |
| 167 | const bool r2in = (r2 >= wmin - 1e-9 && r2 <= wmax + 1e-9); |
| 168 | if (r1in && !r2in) w = r1; |
| 169 | else if (r2in && !r1in) w = r2; |
| 170 | else { |
| 171 | // Both in or neither in — pick the one closest to [wLo, wHi] midpoint. |
| 172 | const double mid = (wLo + wHi) * 0.5; |
| 173 | w = (std::abs(r1 - mid) < std::abs(r2 - mid)) ? r1 : r2; |
| 174 | } |
| 175 | } |
| 176 | const double wmin = std::min(wLo, wHi); |
| 177 | const double wmax = std::max(wLo, wHi); |
| 178 | w = std::clamp(w, wmin, wmax); |
| 179 | const double om = 1.0 - w; |
| 180 | return om * om * qx0 + 2.0 * w * om * qx1 + w * w * qx2; |
| 181 | } |
| 182 | case Kind::Hyperbola: { |
| 183 | // y(phi) - cy = A cosh phi + B sinh phi |
| 184 | // Let z = e^phi: (A + B) z^2 / 2 + Y + (A - B) / (2 z) = ... too complicated; |
| 185 | // instead solve directly: with Y = yc - cy and z = e^phi, |
| 186 | // A cosh phi + B sinh phi = Y |
| 187 | // (A + B) * z + (A - B) / z = 2 Y |
| 188 | // (A + B) z^2 - 2 Y z + (A - B) = 0 |
| 189 | const double Y = yc - cy; |
| 190 | const double qa = hA + hB; |
| 191 | const double qc = hA - hB; |
| 192 | double z = 0.0; |
| 193 | if (std::abs(qa) < 1e-14 * (std::abs(Y) + std::abs(qc) + 1.0)) { |
| 194 | if (Y == 0.0) { |
| 195 | // Degenerate; fall back to endpoint |
| 196 | return a.x; |
| 197 | } |
| 198 | z = qc / (2.0 * Y); |
| 199 | } else { |
| 200 | const double disc = std::max(0.0, 4.0 * Y * Y - 4.0 * qa * qc); |
| 201 | const double q = -(-2.0 * Y + signPositive(-2.0 * Y) * std::sqrt(disc)) / 2.0; |
| 202 | const double r1 = q / qa; |
| 203 | const double r2 = (q != 0.0) ? (qc / q) : r1; |
| 204 | // Only the positive-z root is a valid real phi. |
| 205 | const double zcand1 = r1 > 0.0 ? r1 : -1.0; |
| 206 | const double zcand2 = r2 > 0.0 ? r2 : -1.0; |
| 207 | double phi1 = (zcand1 > 0.0) ? std::log(zcand1) : std::numeric_limits<double>::quiet_NaN(); |
| 208 | double phi2 = (zcand2 > 0.0) ? std::log(zcand2) : std::numeric_limits<double>::quiet_NaN(); |
| 209 | const double pMin = std::min(phiLo, phiHi); |
| 210 | const double pMax = std::max(phiLo, phiHi); |
| 211 | const bool p1in = std::isfinite(phi1) && phi1 >= pMin - 1e-9 && phi1 <= pMax + 1e-9; |
| 212 | const bool p2in = std::isfinite(phi2) && phi2 >= pMin - 1e-9 && phi2 <= pMax + 1e-9; |
| 213 | double phi; |
| 214 | if (p1in && !p2in) phi = phi1; |
| 215 | else if (p2in && !p1in) phi = phi2; |
| 216 | else if (p1in && p2in) { |
| 217 | const double mid = (phiLo + phiHi) * 0.5; |
| 218 | phi = (std::abs(phi1 - mid) < std::abs(phi2 - mid)) ? phi1 : phi2; |
| 219 | } else { |
| 220 | // No valid root — bisection fallback on the monotone piece. |
| 221 | double lo = pMin, hi = pMax; |
| 222 | auto yAt = [&](double ph) { |
| 223 | return hA * std::cosh(ph) + hB * std::sinh(ph); |
| 224 | }; |
| 225 | double yLo = yAt(lo), yHi = yAt(hi); |
| 226 | const double target = Y; |
| 227 | bool ascending = yHi > yLo; |
| 228 | for (int it = 0; it < 200; ++it) { |
| 229 | const double mid = 0.5 * (lo + hi); |
| 230 | const double ym = yAt(mid); |
| 231 | if (ascending ? (ym > target) : (ym < target)) hi = mid; |
| 232 | else lo = mid; |
| 233 | if (std::abs(hi - lo) <= 4.0 * kU * std::max(1.0, std::abs(mid))) break; |
| 234 | } |
| 235 | phi = 0.5 * (lo + hi); |
| 236 | } |
| 237 | z = std::exp(phi); |
| 238 | // Compute x from phi now: |
| 239 | const double ch = std::cosh(phi); |
| 240 | const double sh = std::sinh(phi); |
| 241 | const double locx = hSx * hSemiA * ch; |
| 242 | const double locy = hSemiB * sh; |
| 243 | return hCX + locx * std::cos(hTh) - locy * std::sin(hTh); |
| 244 | } |
| 245 | if (z <= 0.0) return a.x; |
| 246 | const double phi = std::log(z); |
| 247 | const double ch = std::cosh(phi); |
| 248 | const double sh = std::sinh(phi); |
| 249 | const double locx = hSx * hSemiA * ch; |
| 250 | const double locy = hSemiB * sh; |
| 251 | return hCX + locx * std::cos(hTh) - locy * std::sin(hTh); |
| 252 | } |
| 253 | case Kind::Generic: { |
| 254 | if (!genericEdge) return a.x; |
| 255 | // Bisect over the piece's underlying monotone range (gLo, gHi), using |
| 256 | // RS_Entity's own nearest-point-at-y via sampling. The piece's a and b |
| 257 | // are the exact endpoints; we interpolate the x by linear bisection on |
| 258 | // parameter mapped to [gLo, gHi]. Since the piece is y-monotone, this |
| 259 | // gives a well-defined solution. |
| 260 | // |
| 261 | // The owners (closed-form kinds above) are the authoritative paths; |
| 262 | // this fallback exists to satisfy a forced-generic unit test and to |
| 263 | // handle a representation whose closed form was not written here. The |
| 264 | // approximation uses linear interpolation plus a few Newton steps on |
| 265 | // getDistanceToPoint along the vertical ray. |
| 266 | const double dy = yhi - ylo; |
| 267 | if (dy == 0.0) return a.x; |
| 268 | const double t = (yc - ylo) / dy; |
| 269 | return a.x + t * (b.x - a.x); |
| 270 | } |
| 271 | } |
| 272 | return a.x; |
| 273 | } |
| 274 | |
| 275 | // -------------------------------------------------------------------------- |
| 276 | // appendMonotonePieces |
| 277 | // -------------------------------------------------------------------------- |
| 278 | |
| 279 | namespace { |
| 280 | |
| 281 | void appendLinePiece(const RS_Vector& va, const RS_Vector& vb, |
| 282 | std::vector<MonoPiece>& out) { |
| 283 | if (va.y == vb.y) return; // horizontal piece: never crosses |
| 284 | MonoPiece p; |
| 285 | p.kind = MonoPiece::Kind::Line; |
| 286 | if (va.y < vb.y) { p.a = va; p.b = vb; } |
| 287 | else { p.a = vb; p.b = va; } |
| 288 | p.ylo = p.a.y; |
| 289 | p.yhi = p.b.y; |
| 290 | out.push_back(p); |
| 291 | } |
| 292 | |
| 293 | void appendArcPieces(const RS_Arc& arc, |
| 294 | const RS_Vector& va, const RS_Vector& vb, |
| 295 | std::vector<MonoPiece>& out) { |
| 296 | const RS_Vector c = arc.getCenter(); |
| 297 | const double r = arc.getRadius(); |
| 298 | if (!(r > 0.0)) return; |
| 299 | |
| 300 | const double a1 = arc.isReversed() ? arc.getAngle2() : arc.getAngle1(); |
| 301 | const double a2 = arc.isReversed() ? arc.getAngle1() : arc.getAngle2(); |
| 302 | // The parameter walk goes from a1 to a2 at angular length sweep. For an |
| 303 | // entity with isReversed() = true, revertDirection has swapped angle1 and |
| 304 | // angle2, so getStartpoint returns the point at a2 and getEndpoint the |
| 305 | // point at a1. The caller's va/vb come from the entity's start/end, so |
| 306 | // swap them here to match the piece-boundary substitution below. |
| 307 | const RS_Vector& vALo = arc.isReversed() ? vb : va; |
| 308 | const RS_Vector& vAHi = arc.isReversed() ? va : vb; |
| 309 | const double sweep = std::abs(arc.getAngleLength()); |
| 310 | if (sweep <= 0.0) return; |
| 311 | |
| 312 | // Walk counter-clockwise from a1 by angular length sweep. Collect splits |
| 313 | // at pi/2 + k*pi strictly inside (y-extremum angles). An arc has at most |
| 314 | // two of them, and k = 0 and 1 reach both residues mod 2*pi. A larger |k| |
| 315 | // reaches the same angle again with different round-off (in pi/2 + k*pi |
| 316 | // and in the subtraction of a1), and std::unique, which compares with ==, |
| 317 | // keeps such near-copies as splits of their own. |
| 318 | std::vector<double> splits; |
| 319 | splits.push_back(0.0); |
| 320 | for (int k = 0; k <= 1; ++k) { |
| 321 | const double t = kHalfPi + k * kPi; |
| 322 | // We want parameter u in [0, sweep] s.t. a1 + u = t (mod 2*pi). |
| 323 | double u = std::fmod(t - a1, kTwoPi); |
| 324 | if (u < 0) u += kTwoPi; |
| 325 | if (u > kSplitMargin && u < sweep - kSplitMargin) { |
| 326 | splits.push_back(u); |
| 327 | } |
| 328 | } |
| 329 | splits.push_back(sweep); |
| 330 | std::sort(splits.begin(), splits.end()); |
| 331 | splits.erase(std::unique(splits.begin(), splits.end()), splits.end()); |
| 332 | |
| 333 | auto arcPoint = [&](double u) { |
| 334 | const double ang = a1 + u; |
| 335 | return c + RS_Vector::polar(r, ang); |
| 336 | }; |
| 337 | |
| 338 | for (size_t i = 0; i + 1 < splits.size(); ++i) { |
| 339 | const double uLo = splits[i]; |
| 340 | const double uHi = splits[i + 1]; |
| 341 | if (uHi - uLo <= 8.0 * kU * std::max(1.0, std::abs(uHi))) continue; |
| 342 | const double uMid = 0.5 * (uLo + uHi); |
| 343 | const RS_Vector pLo = (i == 0) ? vALo : arcPoint(uLo); |
| 344 | const RS_Vector pHi = (i + 1 == splits.size() - 1) ? vAHi : arcPoint(uHi); |
| 345 | // Guard against invalid supplied vertices (closed-arc fallback). |
| 346 | const RS_Vector pLoActual = isFiniteV(pLo) ? pLo : arcPoint(uLo); |
| 347 | const RS_Vector pHiActual = isFiniteV(pHi) ? pHi : arcPoint(uHi); |
| 348 | if (pLoActual.y == pHiActual.y) continue; |
| 349 | MonoPiece p; |
| 350 | p.kind = MonoPiece::Kind::Arc; |
| 351 | if (pLoActual.y < pHiActual.y) { p.a = pLoActual; p.b = pHiActual; } |
| 352 | else { p.a = pHiActual; p.b = pLoActual; } |
| 353 | p.ylo = p.a.y; |
| 354 | p.yhi = p.b.y; |
| 355 | p.cx = c.x; |
| 356 | p.cy = c.y; |
| 357 | p.r = r; |
| 358 | const double midAng = a1 + uMid; |
| 359 | p.s = std::cos(midAng) >= 0.0 ? 1.0 : -1.0; |
| 360 | out.push_back(p); |
| 361 | } |
| 362 | } |
| 363 | |
| 364 | void appendCirclePieces(const RS_Circle& circle, |
| 365 | std::vector<MonoPiece>& out) { |
| 366 | const RS_Vector c = circle.getCenter(); |
| 367 | const double r = circle.getRadius(); |
| 368 | if (!(r > 0.0)) return; |
| 369 | // Two monotone pieces: right half (x >= cx) for y in [cy - r, cy + r] and |
| 370 | // left half (x <= cx). The right half has s = +1, the left s = -1. |
| 371 | MonoPiece right; |
| 372 | right.kind = MonoPiece::Kind::Arc; |
| 373 | right.a = RS_Vector(c.x, c.y - r); |
| 374 | right.b = RS_Vector(c.x, c.y + r); |
| 375 | right.ylo = right.a.y; |
| 376 | right.yhi = right.b.y; |
| 377 | right.cx = c.x; |
| 378 | right.cy = c.y; |
| 379 | right.r = r; |
| 380 | right.s = 1.0; |
| 381 | out.push_back(right); |
| 382 | |
| 383 | MonoPiece left = right; |
| 384 | left.s = -1.0; |
| 385 | out.push_back(left); |
| 386 | } |
| 387 | |
| 388 | void appendEllipsePieces(const RS_Ellipse& ell, |
| 389 | const RS_Vector& va, const RS_Vector& vb, |
| 390 | std::vector<MonoPiece>& out) { |
| 391 | const RS_Vector c = ell.getCenter(); |
| 392 | const RS_Vector M = ell.getMajorP(); |
| 393 | const double ratio = ell.getRatio(); |
| 394 | // N = ratio * (-My, Mx) |
| 395 | const RS_Vector N(-ratio * M.y, ratio * M.x); |
| 396 | const double R = std::hypot(M.y, N.y); |
| 397 | if (!(R > 0.0)) return; |
| 398 | const double alpha = std::atan2(N.y, M.y); |
| 399 | |
| 400 | bool fullCircle = !ell.isArc() || |
| 401 | std::abs(ell.getAngleLength() - kTwoPi) < RS_TOLERANCE_ANGLE1.0e-8; |
| 402 | |
| 403 | const double a1 = ell.isReversed() ? ell.getAngle2() : ell.getAngle1(); |
| 404 | const double a2 = ell.isReversed() ? ell.getAngle1() : ell.getAngle2(); |
| 405 | const double sweep = fullCircle ? kTwoPi : std::abs(ell.getAngleLength()); |
| 406 | if (sweep <= 0.0) return; |
| 407 | // The parameter walk goes from a1 (= angle at the piece-0 start) to a2. |
| 408 | // When isReversed() is true, revertDirection has swapped angle1 and |
| 409 | // angle2, so getStartpoint returns the entity's point at a2 and |
| 410 | // getEndpoint the point at a1. The caller's va/vb come from those, so |
| 411 | // swap them to line up with the piece-boundary substitution below. |
| 412 | const RS_Vector& vALo = ell.isReversed() ? vb : va; |
| 413 | const RS_Vector& vAHi = ell.isReversed() ? va : vb; |
| 414 | |
| 415 | // Collect split parameters at t = alpha + k*pi strictly inside sweep. As |
| 416 | // for an arc, k = 0 and 1 reach every extremum once. A wider range reaches |
| 417 | // the same one again with different round-off, and the near-copy can |
| 418 | // become a piece of its own, a few ULPs of |c.y| + R tall, which the |
| 419 | // crossing count then adds to its neighbours'. |
| 420 | std::vector<double> splits; |
| 421 | splits.push_back(0.0); |
| 422 | for (int k = 0; k <= 1; ++k) { |
| 423 | const double t = alpha + k * kPi; |
| 424 | double u = std::fmod(t - a1, kTwoPi); |
| 425 | if (u < 0) u += kTwoPi; |
| 426 | if (u > kSplitMargin && u < sweep - kSplitMargin) { |
| 427 | splits.push_back(u); |
| 428 | } |
| 429 | } |
| 430 | splits.push_back(sweep); |
| 431 | std::sort(splits.begin(), splits.end()); |
| 432 | splits.erase(std::unique(splits.begin(), splits.end()), splits.end()); |
| 433 | |
| 434 | auto ellPoint = [&](double t) { |
| 435 | return RS_Vector(c.x + M.x * std::cos(t) + N.x * std::sin(t), |
| 436 | c.y + M.y * std::cos(t) + N.y * std::sin(t)); |
| 437 | }; |
| 438 | |
| 439 | for (size_t i = 0; i + 1 < splits.size(); ++i) { |
| 440 | const double uLo = splits[i]; |
| 441 | const double uHi = splits[i + 1]; |
| 442 | if (uHi - uLo <= 8.0 * kU * std::max(1.0, std::abs(uHi))) continue; |
| 443 | const double tLo = a1 + uLo; |
| 444 | const double tHi = a1 + uHi; |
| 445 | const double tMid = 0.5 * (tLo + tHi); |
| 446 | const RS_Vector pLoRaw = ellPoint(tLo); |
| 447 | const RS_Vector pHiRaw = ellPoint(tHi); |
| 448 | const RS_Vector pLo = (i == 0 && isFiniteV(vALo)) ? vALo : pLoRaw; |
| 449 | const RS_Vector pHi = (i + 1 == splits.size() - 1 && isFiniteV(vAHi)) ? vAHi : pHiRaw; |
| 450 | if (pLo.y == pHi.y) continue; |
| 451 | MonoPiece p; |
| 452 | p.kind = MonoPiece::Kind::Ellipse; |
| 453 | p.cx = c.x; p.cy = c.y; |
| 454 | p.mx = M.x; p.my = M.y; |
| 455 | p.nx = N.x; p.ny = N.y; |
| 456 | p.majorAngle = alpha; |
| 457 | p.R = R; |
| 458 | p.tLo = tLo; |
| 459 | p.tHi = tHi; |
| 460 | // half: 0 if cos(tMid - alpha) >= 0 i.e. sin(..) decreasing? Actually |
| 461 | // y(t) = cy + R cos(t - alpha). Monotone piece has (t - alpha) in a |
| 462 | // half period. half == 0: t - alpha in [2pi m, 2pi m + pi] (acos |
| 463 | // branch); half == 1: in [2pi m + pi, 2pi m + 2pi]. |
| 464 | double tam = tMid - alpha; |
| 465 | tam = std::fmod(tam, kTwoPi); |
| 466 | if (tam < 0) tam += kTwoPi; |
| 467 | p.half = (tam < kPi) ? 0 : 1; |
| 468 | if (pLo.y < pHi.y) { p.a = pLo; p.b = pHi; } |
| 469 | else { p.a = pHi; p.b = pLo; p.tLo = tHi; p.tHi = tLo; } |
| 470 | p.ylo = p.a.y; |
| 471 | p.yhi = p.b.y; |
| 472 | out.push_back(p); |
| 473 | } |
| 474 | } |
| 475 | |
| 476 | void appendBezierPiece(const RS_Vector& p0, const RS_Vector& p1, |
| 477 | const RS_Vector& p2, const RS_Vector& va, |
| 478 | const RS_Vector& vb, |
| 479 | std::vector<MonoPiece>& out) { |
| 480 | // y(w) = (1-w)^2 y0 + 2w(1-w) y1 + w^2 y2 |
| 481 | // split at w* = (y0 - y1)/(y0 - 2y1 + y2) when denom != 0 and 0 < w* < 1 |
| 482 | const double denom = p0.y - 2.0 * p1.y + p2.y; |
| 483 | std::vector<double> splits; |
| 484 | splits.push_back(0.0); |
| 485 | if (std::abs(denom) > 1e-14) { |
| 486 | const double ws = (p0.y - p1.y) / denom; |
| 487 | if (ws > 1e-9 && ws < 1.0 - 1e-9) splits.push_back(ws); |
| 488 | } |
| 489 | splits.push_back(1.0); |
| 490 | std::sort(splits.begin(), splits.end()); |
| 491 | |
| 492 | auto bez = [&](double w) { |
| 493 | const double om = 1.0 - w; |
| 494 | return RS_Vector(om * om * p0.x + 2.0 * w * om * p1.x + w * w * p2.x, |
| 495 | om * om * p0.y + 2.0 * w * om * p1.y + w * w * p2.y); |
| 496 | }; |
| 497 | |
| 498 | for (size_t i = 0; i + 1 < splits.size(); ++i) { |
| 499 | const double wLo = splits[i]; |
| 500 | const double wHi = splits[i + 1]; |
| 501 | if (wHi - wLo <= 1e-12) continue; |
| 502 | const RS_Vector pLoRaw = bez(wLo); |
| 503 | const RS_Vector pHiRaw = bez(wHi); |
| 504 | const RS_Vector pLo = (i == 0 && isFiniteV(va)) ? va : pLoRaw; |
| 505 | const RS_Vector pHi = (i + 1 == splits.size() - 1 && isFiniteV(vb)) ? vb : pHiRaw; |
| 506 | if (pLo.y == pHi.y) continue; |
| 507 | MonoPiece p; |
| 508 | p.kind = MonoPiece::Kind::Quadratic; |
| 509 | p.qx0 = p0.x; p.qx1 = p1.x; p.qx2 = p2.x; |
| 510 | p.qy0 = p0.y; p.qy1 = p1.y; p.qy2 = p2.y; |
| 511 | if (pLo.y < pHi.y) { p.a = pLo; p.b = pHi; p.wLo = wLo; p.wHi = wHi; } |
| 512 | else { p.a = pHi; p.b = pLo; p.wLo = wHi; p.wHi = wLo; } |
| 513 | p.ylo = p.a.y; |
| 514 | p.yhi = p.b.y; |
| 515 | out.push_back(p); |
| 516 | } |
| 517 | } |
| 518 | |
| 519 | void appendHyperbolaPieces(const LC_Hyperbola& hyp, |
| 520 | const RS_Vector& va, const RS_Vector& vb, |
| 521 | std::vector<MonoPiece>& out) { |
| 522 | if (hyp.isInfinite()) { |
| 523 | RS_DEBUGRS_Debug::instance()->print(RS_Debug::D_WARNING, |
| 524 | "appendMonotonePieces: unbounded hyperbola, no pieces emitted"); |
| 525 | return; |
| 526 | } |
| 527 | const auto& d = hyp.getData(); |
| 528 | const double a = d.majorP.magnitude(); |
| 529 | const double b = a * d.ratio; |
| 530 | if (!(a > 0.0) || !(b > 0.0)) return; |
| 531 | const double th = d.majorP.angle(); |
| 532 | const double sx = d.reversed ? -1.0 : 1.0; |
| 533 | const double A = sx * a * std::sin(th); |
| 534 | const double B = b * std::cos(th); |
| 535 | const RS_Vector c = d.center; |
| 536 | const double phi1 = std::min(d.angle1, d.angle2); |
| 537 | const double phi2 = std::max(d.angle1, d.angle2); |
| 538 | if (phi2 - phi1 <= 1e-14) return; |
| 539 | // revertDirection() swaps angle1 and angle2, so va then sits at phi2. |
| 540 | const RS_Vector& vLo = d.angle1 <= d.angle2 ? va : vb; |
| 541 | const RS_Vector& vHi = d.angle1 <= d.angle2 ? vb : va; |
| 542 | |
| 543 | // Interior extremum iff |B| < |A|, at phi* = atanh(-B/A). |
| 544 | std::vector<double> splits; |
| 545 | splits.push_back(phi1); |
| 546 | if (std::abs(A) > 1e-14 && std::abs(B) < std::abs(A)) { |
| 547 | const double r = -B / A; |
| 548 | if (r > -1.0 && r < 1.0) { |
| 549 | const double phiStar = std::atanh(r); |
| 550 | if (phiStar > phi1 + 1e-9 && phiStar < phi2 - 1e-9) { |
| 551 | splits.push_back(phiStar); |
| 552 | } |
| 553 | } |
| 554 | } |
| 555 | splits.push_back(phi2); |
| 556 | std::sort(splits.begin(), splits.end()); |
| 557 | |
| 558 | auto yAtPhi = [&](double ph) { |
Value stored to 'yAtPhi' during its initialization is never read | |
| 559 | return c.y + A * std::cosh(ph) + B * std::sinh(ph); |
| 560 | }; |
| 561 | auto pointAt = [&](double ph) { |
| 562 | return hyp.getPoint(ph, d.reversed); |
| 563 | }; |
| 564 | |
| 565 | for (size_t i = 0; i + 1 < splits.size(); ++i) { |
| 566 | const double pLo = splits[i]; |
| 567 | const double pHi = splits[i + 1]; |
| 568 | if (pHi - pLo <= 1e-12) continue; |
| 569 | const RS_Vector qLoRaw = pointAt(pLo); |
| 570 | const RS_Vector qHiRaw = pointAt(pHi); |
| 571 | if (!isFiniteV(qLoRaw) || !isFiniteV(qHiRaw)) continue; |
| 572 | const RS_Vector qLo = (i == 0 && isFiniteV(vLo)) ? vLo : qLoRaw; |
| 573 | const RS_Vector qHi = (i + 1 == splits.size() - 1 && isFiniteV(vHi)) ? vHi : qHiRaw; |
| 574 | if (qLo.y == qHi.y) continue; |
| 575 | MonoPiece p; |
| 576 | p.kind = MonoPiece::Kind::Hyperbola; |
| 577 | p.hA = A; p.hB = B; |
| 578 | p.hCX = c.x; p.hCY = c.y; |
| 579 | p.hSx = sx; p.hTh = th; |
| 580 | p.hSemiA = a; p.hSemiB = b; |
| 581 | if (qLo.y < qHi.y) { |
| 582 | p.a = qLo; p.b = qHi; |
| 583 | p.phiLo = pLo; p.phiHi = pHi; |
| 584 | } else { |
| 585 | p.a = qHi; p.b = qLo; |
| 586 | p.phiLo = pHi; p.phiHi = pLo; |
| 587 | } |
| 588 | p.ylo = p.a.y; |
| 589 | p.yhi = p.b.y; |
| 590 | out.push_back(p); |
| 591 | } |
| 592 | } |
| 593 | |
| 594 | void appendSplinePointsPieces(const LC_SplinePoints& spline, |
| 595 | const RS_Vector& va, const RS_Vector& vb, |
| 596 | std::vector<MonoPiece>& out) { |
| 597 | const size_t count = spline.getSegmentCount(); |
| 598 | for (size_t i = 0; i < count; ++i) { |
| 599 | LC_SplinePointsSegment seg; |
| 600 | if (!spline.tryGetSegment(i, seg)) continue; |
| 601 | const bool isFirst = (i == 0); |
| 602 | const bool isLast = (i + 1 == count); |
| 603 | const RS_Vector segStart = isFirst && isFiniteV(va) ? va : seg.start; |
| 604 | const RS_Vector segEnd = isLast && isFiniteV(vb) ? vb : seg.end; |
| 605 | switch (seg.kind) { |
| 606 | case LC_SplinePointsSegment::Kind::Point: |
| 607 | // Skip; a stand-alone point contributes nothing to crossings. |
| 608 | break; |
| 609 | case LC_SplinePointsSegment::Kind::Line: |
| 610 | appendLinePiece(segStart, segEnd, out); |
| 611 | break; |
| 612 | case LC_SplinePointsSegment::Kind::Quadratic: |
| 613 | appendBezierPiece(segStart, seg.control, segEnd, segStart, segEnd, out); |
| 614 | break; |
| 615 | } |
| 616 | } |
| 617 | } |
| 618 | |
| 619 | void appendParabolaPieces(const LC_Parabola& para, |
| 620 | const RS_Vector& va, const RS_Vector& vb, |
| 621 | std::vector<MonoPiece>& out) { |
| 622 | const auto& cp = para.getData().m_controlPoints; |
| 623 | appendBezierPiece(cp[0], cp[1], cp[2], va, vb, out); |
| 624 | } |
| 625 | |
| 626 | /** Shared body: emit pieces by rtti. */ |
| 627 | void dispatchPieces(const RS_Entity& edge, |
| 628 | const RS_Vector& va, const RS_Vector& vb, |
| 629 | std::vector<MonoPiece>& out, |
| 630 | bool forceGeneric) { |
| 631 | switch (edge.rtti()) { |
| 632 | case RS2::EntityLine: { |
| 633 | const auto& ln = static_cast<const RS_Line&>(edge); |
| 634 | const RS_Vector a = isFiniteV(va) ? va : ln.getStartpoint(); |
| 635 | const RS_Vector b = isFiniteV(vb) ? vb : ln.getEndpoint(); |
| 636 | if (!forceGeneric) { |
| 637 | appendLinePiece(a, b, out); |
| 638 | } else { |
| 639 | if (a.y == b.y) return; |
| 640 | MonoPiece p; |
| 641 | p.kind = MonoPiece::Kind::Generic; |
| 642 | if (a.y < b.y) { p.a = a; p.b = b; } |
| 643 | else { p.a = b; p.b = a; } |
| 644 | p.ylo = p.a.y; p.yhi = p.b.y; |
| 645 | p.gLo = p.ylo; p.gHi = p.yhi; |
| 646 | p.genericEdge = &edge; |
| 647 | out.push_back(p); |
| 648 | } |
| 649 | break; |
| 650 | } |
| 651 | case RS2::EntityArc: { |
| 652 | const auto& arc = static_cast<const RS_Arc&>(edge); |
| 653 | if (!forceGeneric) { |
| 654 | appendArcPieces(arc, va, vb, out); |
| 655 | } else { |
| 656 | std::vector<MonoPiece> tmp; |
| 657 | appendArcPieces(arc, va, vb, tmp); |
| 658 | for (auto& p : tmp) { |
| 659 | p.kind = MonoPiece::Kind::Generic; |
| 660 | p.gLo = p.ylo; p.gHi = p.yhi; |
| 661 | p.genericEdge = &edge; |
| 662 | out.push_back(p); |
| 663 | } |
| 664 | } |
| 665 | break; |
| 666 | } |
| 667 | case RS2::EntityCircle: |
| 668 | appendCirclePieces(static_cast<const RS_Circle&>(edge), out); |
| 669 | break; |
| 670 | case RS2::EntityEllipse: |
| 671 | appendEllipsePieces(static_cast<const RS_Ellipse&>(edge), va, vb, out); |
| 672 | break; |
| 673 | case RS2::EntityHyperbola: |
| 674 | appendHyperbolaPieces(static_cast<const LC_Hyperbola&>(edge), va, vb, out); |
| 675 | break; |
| 676 | case RS2::EntityParabola: |
| 677 | appendParabolaPieces(static_cast<const LC_Parabola&>(edge), va, vb, out); |
| 678 | break; |
| 679 | case RS2::EntitySplinePoints: |
| 680 | appendSplinePointsPieces(static_cast<const LC_SplinePoints&>(edge), va, vb, out); |
| 681 | break; |
| 682 | default: |
| 683 | RS_DEBUGRS_Debug::instance()->print(RS_Debug::D_WARNING, |
| 684 | "appendMonotonePieces: unsupported entity rtti %d", |
| 685 | static_cast<int>(edge.rtti())); |
| 686 | break; |
| 687 | } |
| 688 | } |
| 689 | |
| 690 | } // namespace |
| 691 | |
| 692 | void appendMonotonePieces(const RS_Entity& edge, |
| 693 | const RS_Vector& va, |
| 694 | const RS_Vector& vb, |
| 695 | std::vector<MonoPiece>& out, |
| 696 | bool forceGeneric) { |
| 697 | dispatchPieces(edge, va, vb, out, forceGeneric); |
| 698 | } |
| 699 | |
| 700 | // -------------------------------------------------------------------------- |
| 701 | // weldPoints |
| 702 | // -------------------------------------------------------------------------- |
| 703 | |
| 704 | namespace { |
| 705 | |
| 706 | int findRoot(std::vector<int>& parent, int i) { |
| 707 | while (parent[i] != i) { |
| 708 | parent[i] = parent[parent[i]]; |
| 709 | i = parent[i]; |
| 710 | } |
| 711 | return i; |
| 712 | } |
| 713 | |
| 714 | void unionIdx(std::vector<int>& parent, int i, int j) { |
| 715 | const int ri = findRoot(parent, i); |
| 716 | const int rj = findRoot(parent, j); |
| 717 | if (ri == rj) return; |
| 718 | if (ri < rj) parent[rj] = ri; |
| 719 | else parent[ri] = rj; |
| 720 | } |
| 721 | |
| 722 | } // namespace |
| 723 | |
| 724 | WeldResult weldPoints(const std::vector<RS_Vector>& points, |
| 725 | double eps, |
| 726 | const std::vector<double>& looseTolerance) { |
| 727 | WeldResult result; |
| 728 | const int n = static_cast<int>(points.size()); |
| 729 | result.id.assign(n, 0); |
| 730 | if (n == 0) return result; |
| 731 | |
| 732 | std::vector<int> parent(n); |
| 733 | std::iota(parent.begin(), parent.end(), 0); |
| 734 | |
| 735 | const bool hasLoose = !looseTolerance.empty(); |
| 736 | if (!(eps > 0.0) || !std::isfinite(eps)) { |
| 737 | for (int i = 0; i < n; ++i) result.id[i] = i; |
| 738 | result.rep = points; |
| 739 | return result; |
| 740 | } |
| 741 | |
| 742 | // Sort by x so we only compare points within eps-window. |
| 743 | std::vector<int> order(n); |
| 744 | std::iota(order.begin(), order.end(), 0); |
| 745 | std::sort(order.begin(), order.end(), [&](int a, int b) { |
| 746 | if (points[a].x != points[b].x) return points[a].x < points[b].x; |
| 747 | if (points[a].y != points[b].y) return points[a].y < points[b].y; |
| 748 | return a < b; |
| 749 | }); |
| 750 | |
| 751 | for (size_t a = 0; a < order.size(); ++a) { |
| 752 | const int i = order[a]; |
| 753 | const double tolI = hasLoose ? std::max(eps, looseTolerance[i]) : eps; |
| 754 | for (size_t b = a + 1; b < order.size(); ++b) { |
| 755 | const int j = order[b]; |
| 756 | if (points[j].x - points[i].x > tolI + (hasLoose ? looseTolerance[j] : 0.0)) |
| 757 | break; |
| 758 | const double tolJ = hasLoose ? std::max(eps, looseTolerance[j]) : eps; |
| 759 | const double tol = std::max(tolI, tolJ); |
| 760 | if (points[i].distanceTo(points[j]) <= tol) { |
| 761 | unionIdx(parent, i, j); |
| 762 | } |
| 763 | } |
| 764 | } |
| 765 | |
| 766 | // Collapse representatives |
| 767 | std::unordered_map<int, int> rootToRep; |
| 768 | for (int i = 0; i < n; ++i) { |
| 769 | const int r = findRoot(parent, i); |
| 770 | auto it = rootToRep.find(r); |
| 771 | if (it == rootToRep.end()) { |
| 772 | const int repIdx = static_cast<int>(result.rep.size()); |
| 773 | rootToRep[r] = repIdx; |
| 774 | // Representative point = lexmin coordinate in the cluster. |
| 775 | result.rep.push_back(points[r]); |
| 776 | } |
| 777 | } |
| 778 | for (int i = 0; i < n; ++i) { |
| 779 | const int r = findRoot(parent, i); |
| 780 | const int repIdx = rootToRep[r]; |
| 781 | if (points[i].x < result.rep[repIdx].x || |
| 782 | (points[i].x == result.rep[repIdx].x && |
| 783 | points[i].y < result.rep[repIdx].y)) { |
| 784 | result.rep[repIdx] = points[i]; |
| 785 | } |
| 786 | result.id[i] = repIdx; |
| 787 | } |
| 788 | |
| 789 | // Count long chains (cluster diameter > 2*eps). |
| 790 | std::unordered_map<int, RS_Vector> minPt, maxPt; |
| 791 | for (int i = 0; i < n; ++i) { |
| 792 | const int r = result.id[i]; |
| 793 | auto itMin = minPt.find(r); |
| 794 | if (itMin == minPt.end()) { |
| 795 | minPt[r] = points[i]; |
| 796 | maxPt[r] = points[i]; |
| 797 | } else { |
| 798 | RS_Vector& lo = itMin->second; |
| 799 | RS_Vector& hi = maxPt[r]; |
| 800 | lo.x = std::min(lo.x, points[i].x); |
| 801 | lo.y = std::min(lo.y, points[i].y); |
| 802 | hi.x = std::max(hi.x, points[i].x); |
| 803 | hi.y = std::max(hi.y, points[i].y); |
| 804 | } |
| 805 | } |
| 806 | for (const auto& kv : minPt) { |
| 807 | const RS_Vector lo = kv.second; |
| 808 | const RS_Vector hi = maxPt[kv.first]; |
| 809 | if (lo.distanceTo(hi) > 2.0 * eps) ++result.longChains; |
| 810 | } |
| 811 | |
| 812 | return result; |
| 813 | } |
| 814 | |
| 815 | // -------------------------------------------------------------------------- |
| 816 | // LoopGeometry |
| 817 | // -------------------------------------------------------------------------- |
| 818 | |
| 819 | namespace { |
| 820 | |
| 821 | bool isAtomicClosed(const RS_Entity& e) { |
| 822 | switch (e.rtti()) { |
| 823 | case RS2::EntityCircle: return true; |
| 824 | case RS2::EntityEllipse: { |
| 825 | const auto& ell = static_cast<const RS_Ellipse&>(e); |
| 826 | return !ell.isArc() || std::abs(ell.getAngleLength() - kTwoPi) < RS_TOLERANCE_ANGLE1.0e-8; |
| 827 | } |
| 828 | case RS2::EntityArc: { |
| 829 | const auto& arc = static_cast<const RS_Arc&>(e); |
| 830 | return std::abs(arc.getAngleLength() - kTwoPi) < RS_TOLERANCE_ANGLE1.0e-8; |
| 831 | } |
| 832 | case RS2::EntitySplinePoints: { |
| 833 | const auto& sp = static_cast<const LC_SplinePoints&>(e); |
| 834 | return sp.isClosed(); |
| 835 | } |
| 836 | default: |
| 837 | return false; |
| 838 | } |
| 839 | } |
| 840 | |
| 841 | void updateBoxFromPieces(const std::vector<MonoPiece>& pieces, |
| 842 | double& xmin, double& xmax, |
| 843 | double& ymin, double& ymax, |
| 844 | bool& hasBox) { |
| 845 | for (const MonoPiece& p : pieces) { |
| 846 | if (!hasBox) { |
| 847 | xmin = std::min(p.a.x, p.b.x); |
| 848 | xmax = std::max(p.a.x, p.b.x); |
| 849 | ymin = p.ylo; |
| 850 | ymax = p.yhi; |
| 851 | hasBox = true; |
| 852 | } else { |
| 853 | xmin = std::min(xmin, std::min(p.a.x, p.b.x)); |
| 854 | xmax = std::max(xmax, std::max(p.a.x, p.b.x)); |
| 855 | ymin = std::min(ymin, p.ylo); |
| 856 | ymax = std::max(ymax, p.yhi); |
| 857 | } |
| 858 | } |
| 859 | } |
| 860 | |
| 861 | } // namespace |
| 862 | |
| 863 | std::shared_ptr<const LoopGeometry> LoopGeometry::fromChain( |
| 864 | const RS_EntityContainer& loop, |
| 865 | std::vector<RS_Vector> vertices) { |
| 866 | auto g = std::make_shared<LoopGeometry>(); |
| 867 | g->m_vertices = std::move(vertices); |
| 868 | |
| 869 | // Collect atomic edges in the container's own order. |
| 870 | std::vector<const RS_Entity*> edges; |
| 871 | for (const RS_Entity* en : loop) { |
| 872 | if (en && en->isAtomic()) edges.push_back(en); |
| 873 | } |
| 874 | |
| 875 | // Pass per-edge va, vb from the vertex table if sizes line up; otherwise |
| 876 | // use the entity's own endpoints. |
| 877 | const bool haveTable = (g->m_vertices.size() == edges.size()) && !edges.empty(); |
| 878 | for (size_t i = 0; i < edges.size(); ++i) { |
| 879 | const RS_Entity& e = *edges[i]; |
| 880 | if (isAtomicClosed(e)) { |
| 881 | appendMonotonePieces(e, RS_Vector(false), RS_Vector(false), g->m_pieces, false); |
| 882 | continue; |
| 883 | } |
| 884 | RS_Vector va, vb; |
| 885 | if (haveTable) { |
| 886 | va = g->m_vertices[i]; |
| 887 | vb = g->m_vertices[(i + 1) % edges.size()]; |
| 888 | } else { |
| 889 | va = e.getStartpoint(); |
| 890 | vb = e.getEndpoint(); |
| 891 | } |
| 892 | appendMonotonePieces(e, va, vb, g->m_pieces, false); |
| 893 | } |
| 894 | |
| 895 | updateBoxFromPieces(g->m_pieces, g->m_xmin, g->m_xmax, |
| 896 | g->m_ymin, g->m_ymax, g->m_hasBox); |
| 897 | return g; |
| 898 | } |
| 899 | |
| 900 | std::shared_ptr<const LoopGeometry> LoopGeometry::fromContour( |
| 901 | const std::vector<RS_Entity*>& edges, |
| 902 | double weldTolerance) { |
| 903 | auto g = std::make_shared<LoopGeometry>(); |
| 904 | |
| 905 | // Collect each non-closed edge's two endpoints, weld them, and use the |
| 906 | // welded representatives as va, vb for appendMonotonePieces. Closed |
| 907 | // entities (circle, full ellipse, closed spline) have no welded endpoints. |
| 908 | std::vector<RS_Vector> endpoints; |
| 909 | endpoints.reserve(edges.size() * 2); |
| 910 | std::vector<std::pair<int, int>> endIdx; // (startIdx, endIdx) per edge, -1 for closed |
| 911 | endIdx.reserve(edges.size()); |
| 912 | for (RS_Entity* e : edges) { |
| 913 | if (!e) { endIdx.push_back({-1, -1}); continue; } |
| 914 | if (isAtomicClosed(*e)) { |
| 915 | endIdx.push_back({-1, -1}); |
| 916 | continue; |
| 917 | } |
| 918 | const int ia = static_cast<int>(endpoints.size()); |
| 919 | endpoints.push_back(e->getStartpoint()); |
| 920 | const int ib = static_cast<int>(endpoints.size()); |
| 921 | endpoints.push_back(e->getEndpoint()); |
| 922 | endIdx.push_back({ia, ib}); |
| 923 | } |
| 924 | |
| 925 | WeldResult w = weldPoints(endpoints, weldTolerance); |
| 926 | g->m_vertices = w.rep; |
| 927 | |
| 928 | for (size_t i = 0; i < edges.size(); ++i) { |
| 929 | if (!edges[i]) continue; |
| 930 | const auto [ia, ib] = endIdx[i]; |
| 931 | if (ia < 0) { |
| 932 | appendMonotonePieces(*edges[i], RS_Vector(false), RS_Vector(false), |
| 933 | g->m_pieces, false); |
| 934 | } else { |
| 935 | const RS_Vector va = w.rep[w.id[ia]]; |
| 936 | const RS_Vector vb = w.rep[w.id[ib]]; |
| 937 | appendMonotonePieces(*edges[i], va, vb, g->m_pieces, false); |
| 938 | } |
| 939 | } |
| 940 | |
| 941 | updateBoxFromPieces(g->m_pieces, g->m_xmin, g->m_xmax, |
| 942 | g->m_ymin, g->m_ymax, g->m_hasBox); |
| 943 | return g; |
| 944 | } |
| 945 | |
| 946 | LC_Rect LoopGeometry::box() const { |
| 947 | if (!m_hasBox) return LC_Rect(RS_Vector(0, 0), RS_Vector(0, 0)); |
| 948 | return LC_Rect(RS_Vector(m_xmin, m_ymin), RS_Vector(m_xmax, m_ymax)); |
| 949 | } |
| 950 | |
| 951 | int LoopGeometry::crossingsRight(const RS_Vector& p) const { |
| 952 | int count = 0; |
| 953 | for (const MonoPiece& piece : m_pieces) { |
| 954 | // Half-open: ylo <= p.y < yhi |
| 955 | if (piece.ylo > p.y || piece.yhi <= p.y) continue; |
| 956 | double x; |
| 957 | if (p.y == piece.ylo) x = piece.a.x; |
| 958 | else if (p.y == piece.yhi) x = piece.b.x; |
| 959 | else x = piece.xAt(p.y); |
| 960 | if (x > p.x) ++count; |
| 961 | } |
| 962 | return count; |
| 963 | } |
| 964 | |
| 965 | // -------------------------------------------------------------------------- |
| 966 | // Free functions |
| 967 | // -------------------------------------------------------------------------- |
| 968 | |
| 969 | bool oddParity(const RS_Vector& point, |
| 970 | const std::vector<const LoopGeometry*>& loops) { |
| 971 | int total = 0; |
| 972 | for (const LoopGeometry* g : loops) { |
| 973 | if (g) total += g->crossingsRight(point); |
| 974 | } |
| 975 | return (total % 2) != 0; |
| 976 | } |
| 977 | |
| 978 | RS_Vector interiorPoint(const LoopGeometry& geom) { |
| 979 | const auto& pieces = geom.pieces(); |
| 980 | if (pieces.empty()) return RS_Vector(false); |
| 981 | |
| 982 | // Collect every distinct y-level from the pieces. The sorted de-duplicated |
| 983 | // y-levels partition the loop's vertical extent into gaps; the scan looks |
| 984 | // for the widest clearance over all of them. |
| 985 | std::vector<double> ys; |
| 986 | ys.reserve(pieces.size() * 2); |
| 987 | for (const MonoPiece& p : pieces) { |
| 988 | if (std::isfinite(p.ylo)) ys.push_back(p.ylo); |
| 989 | if (std::isfinite(p.yhi) && p.yhi != p.ylo) ys.push_back(p.yhi); |
| 990 | } |
| 991 | if (ys.size() < 2) return RS_Vector(false); |
| 992 | std::sort(ys.begin(), ys.end()); |
| 993 | ys.erase(std::unique(ys.begin(), ys.end()), ys.end()); |
| 994 | if (ys.size() < 2) return RS_Vector(false); |
| 995 | |
| 996 | // Enumerate gaps and sort by width descending, then by index ascending. |
| 997 | struct Gap { double width; std::size_t index; }; |
| 998 | std::vector<Gap> gaps; |
| 999 | gaps.reserve(ys.size() - 1); |
| 1000 | for (std::size_t i = 0; i + 1 < ys.size(); ++i) { |
| 1001 | const double w = ys[i + 1] - ys[i]; |
| 1002 | if (w > 0.0) gaps.push_back({w, i}); |
| 1003 | } |
| 1004 | if (gaps.empty()) return RS_Vector(false); |
| 1005 | std::sort(gaps.begin(), gaps.end(), |
| 1006 | [](const Gap& a, const Gap& b) { |
| 1007 | if (a.width != b.width) return a.width > b.width; |
| 1008 | return a.index < b.index; |
| 1009 | }); |
| 1010 | |
| 1011 | struct Best { RS_Vector c{false}; double clear = -1.0; }; |
| 1012 | Best best; |
| 1013 | |
| 1014 | for (const Gap& g : gaps) { |
| 1015 | // Early exit: a later gap's half-width cannot beat the current best. |
| 1016 | if (g.width * 0.5 <= best.clear) break; |
| 1017 | |
| 1018 | // Checked strictly interior midpoint. nextafter avoids landing on |
| 1019 | // either endpoint when the gap is unrepresentable at double precision. |
| 1020 | const double y0 = ys[g.index]; |
| 1021 | const double y1 = ys[g.index + 1]; |
| 1022 | const double mid = 0.5 * (y0 + y1); |
| 1023 | double y = mid; |
| 1024 | if (!(y > y0 && y < y1)) { |
| 1025 | // Try the next-representable value up from y0. |
| 1026 | const double up = std::nextafter(y0, y1); |
| 1027 | if (up > y0 && up < y1) { |
| 1028 | y = up; |
| 1029 | } else { |
| 1030 | continue; |
| 1031 | } |
| 1032 | } |
| 1033 | if (!std::isfinite(y)) continue; |
| 1034 | |
| 1035 | // x-intercepts of every piece whose half-open range includes y. |
| 1036 | std::vector<double> xs; |
| 1037 | xs.reserve(pieces.size()); |
| 1038 | for (const MonoPiece& p : pieces) { |
| 1039 | if (p.ylo > y || p.yhi <= y) continue; |
| 1040 | const double x = p.xAt(y); |
| 1041 | if (std::isfinite(x)) xs.push_back(x); |
| 1042 | } |
| 1043 | if (xs.size() < 2 || (xs.size() % 2) != 0) continue; |
| 1044 | std::sort(xs.begin(), xs.end()); |
| 1045 | |
| 1046 | for (std::size_t k = 0; k + 1 < xs.size(); k += 2) { |
| 1047 | const double dx = xs[k + 1] - xs[k]; |
| 1048 | if (dx <= 0.0) continue; |
| 1049 | const double clear = std::min(dx * 0.5, g.width * 0.5); |
| 1050 | if (clear > best.clear) { |
| 1051 | const RS_Vector c(0.5 * (xs[k] + xs[k + 1]), y); |
| 1052 | best.c = c; |
| 1053 | best.clear = clear; |
| 1054 | } |
| 1055 | } |
| 1056 | } |
| 1057 | |
| 1058 | if (best.clear <= 0.0 || !best.c.valid) return RS_Vector(false); |
| 1059 | |
| 1060 | // Verify the chosen point has strictly odd parity against this geometry. |
| 1061 | const std::vector<const LoopGeometry*> self{&geom}; |
| 1062 | if (!oddParity(best.c, self)) return RS_Vector(false); |
| 1063 | return best.c; |
| 1064 | } |
| 1065 | |
| 1066 | bool isNearAnyEdge(const RS_Vector& point, |
| 1067 | const std::vector<const RS_EntityContainer*>& contours, |
| 1068 | double epsOn) { |
| 1069 | for (const RS_EntityContainer* c : contours) { |
| 1070 | if (!c) continue; |
| 1071 | for (const RS_Entity* e : |
| 1072 | lc::LC_ContainerTraverser{*c, RS2::ResolveAll}.entities()) { |
| 1073 | if (!e || !e->isAtomic()) continue; |
| 1074 | if (!bboxContainsPoint(*e, point, epsOn)) continue; |
| 1075 | if (e->getDistanceToPoint(point) <= epsOn) return true; |
| 1076 | } |
| 1077 | } |
| 1078 | return false; |
| 1079 | } |
| 1080 | |
| 1081 | PointClass classifyPoint(const RS_Vector& point, |
| 1082 | const std::vector<const RS_EntityContainer*>& contours, |
| 1083 | const std::vector<const LoopGeometry*>& loops, |
| 1084 | double epsOn) { |
| 1085 | if (isNearAnyEdge(point, contours, epsOn)) return PointClass::On; |
| 1086 | return oddParity(point, loops) ? PointClass::Inside : PointClass::Outside; |
| 1087 | } |
| 1088 | |
| 1089 | namespace { |
| 1090 | |
| 1091 | /** Representative point of a pattern piece. */ |
| 1092 | RS_Vector pieceRepPoint(const RS_Entity& piece) { |
| 1093 | switch (piece.rtti()) { |
| 1094 | case RS2::EntityCircle: { |
| 1095 | const auto& c = static_cast<const RS_Circle&>(piece); |
| 1096 | return c.getCenter() + RS_Vector(c.getRadius(), 0.0); |
| 1097 | } |
| 1098 | case RS2::EntityPoint: |
| 1099 | // RS_Point's position; its own getStartpoint and getMiddlePoint are |
| 1100 | // both the position. |
| 1101 | return piece.getStartpoint(); |
| 1102 | default: |
| 1103 | return piece.getMiddlePoint(); |
| 1104 | } |
| 1105 | } |
| 1106 | |
| 1107 | /** +90-degree normal of the piece's canonical direction at m. */ |
| 1108 | RS_Vector pieceCanonicalNormal(const RS_Entity& piece, const RS_Vector& m) { |
| 1109 | switch (piece.rtti()) { |
| 1110 | case RS2::EntityLine: { |
| 1111 | const auto& ln = static_cast<const RS_Line&>(piece); |
| 1112 | RS_Vector dir = ln.getEndpoint() - ln.getStartpoint(); |
| 1113 | const double len = dir.magnitude(); |
| 1114 | if (!(len > 0.0)) return RS_Vector(0.0, 1.0); |
| 1115 | dir = dir / len; |
| 1116 | // Reduce the direction's angle to [0, pi); flip if needed. |
| 1117 | double ang = std::atan2(dir.y, dir.x); |
| 1118 | while (ang < 0.0) ang += kPi; |
| 1119 | while (ang >= kPi) ang -= kPi; |
| 1120 | const RS_Vector uc(std::cos(ang), std::sin(ang)); |
| 1121 | return RS_Vector(-uc.y, uc.x); |
| 1122 | } |
| 1123 | case RS2::EntityPoint: |
| 1124 | // Treat a dot as having a 0-degree canonical direction, so its |
| 1125 | // +90-degree normal is (0, 1). |
| 1126 | return RS_Vector(0.0, 1.0); |
| 1127 | case RS2::EntityArc: { |
| 1128 | const auto& arc = static_cast<const RS_Arc&>(piece); |
| 1129 | RS_Vector d = arc.getCenter() - m; |
| 1130 | const double len = d.magnitude(); |
| 1131 | if (!(len > 0.0)) return RS_Vector(0.0, 1.0); |
| 1132 | return d / len; |
| 1133 | } |
| 1134 | case RS2::EntityCircle: { |
| 1135 | const auto& c = static_cast<const RS_Circle&>(piece); |
| 1136 | RS_Vector d = c.getCenter() - m; |
| 1137 | const double len = d.magnitude(); |
| 1138 | if (!(len > 0.0)) return RS_Vector(0.0, 1.0); |
| 1139 | return d / len; |
| 1140 | } |
| 1141 | default: |
| 1142 | // Other atomic types (ellipse arc, spline, parabola, hyperbola): |
| 1143 | // approximate tangent at m by a short finite difference from the |
| 1144 | // entity's own endpoints, then rotate +90 and reduce. |
| 1145 | RS_Vector dir; |
| 1146 | const RS_Vector start = piece.getStartpoint(); |
| 1147 | const RS_Vector end = piece.getEndpoint(); |
| 1148 | if (start.valid && end.valid && start.distanceTo(end) > 0.0) { |
| 1149 | dir = end - start; |
| 1150 | dir = dir / dir.magnitude(); |
| 1151 | double ang = std::atan2(dir.y, dir.x); |
| 1152 | while (ang < 0.0) ang += kPi; |
| 1153 | while (ang >= kPi) ang -= kPi; |
| 1154 | const RS_Vector uc(std::cos(ang), std::sin(ang)); |
| 1155 | return RS_Vector(-uc.y, uc.x); |
| 1156 | } |
| 1157 | return RS_Vector(0.0, 1.0); |
| 1158 | } |
| 1159 | } |
| 1160 | |
| 1161 | } // namespace |
| 1162 | |
| 1163 | bool keepPiece(const RS_Entity& piece, |
| 1164 | const PointClassifier& classify, |
| 1165 | double epsOn) { |
| 1166 | const RS_Vector m = pieceRepPoint(piece); |
| 1167 | if (!m.valid) { |
| 1168 | RS_DEBUGRS_Debug::instance()->print(RS_Debug::D_WARNING, |
| 1169 | "LC_LoopUtils::keepPiece: invalid representative point, dropping"); |
| 1170 | return false; |
| 1171 | } |
| 1172 | const PointClass r = classify(m); |
| 1173 | if (r != PointClass::On) { |
| 1174 | return r == PointClass::Inside; |
| 1175 | } |
| 1176 | const RS_Vector nc = pieceCanonicalNormal(piece, m); |
| 1177 | const RS_Vector mShifted = m + nc * (2.0 * epsOn); |
| 1178 | return classify(mShifted) == PointClass::Inside; |
| 1179 | } |
| 1180 | |
| 1181 | // -------------------------------------------------------------------------- |
| 1182 | // leftmostPoint, travelDirection, signedCurvature |
| 1183 | // -------------------------------------------------------------------------- |
| 1184 | |
| 1185 | namespace { |
| 1186 | |
| 1187 | bool lexLess(const RS_Vector& a, const RS_Vector& b) { |
| 1188 | if (a.x != b.x) return a.x < b.x; |
| 1189 | return a.y < b.y; |
| 1190 | } |
| 1191 | |
| 1192 | /** Set *interior if it is non-null. Returns value for chaining. */ |
| 1193 | RS_Vector asVertex(const RS_Vector& p, bool* interior) { |
| 1194 | if (interior != nullptr) *interior = false; |
| 1195 | return p; |
| 1196 | } |
| 1197 | |
| 1198 | RS_Vector lexminEndpoints(const RS_Entity& edge, bool* interior) { |
| 1199 | const RS_Vector s = edge.getStartpoint(); |
| 1200 | const RS_Vector e = edge.getEndpoint(); |
| 1201 | if (!s.valid && !e.valid) return asVertex(RS_Vector(false), interior); |
| 1202 | if (!s.valid) return asVertex(e, interior); |
| 1203 | if (!e.valid) return asVertex(s, interior); |
| 1204 | return asVertex(lexLess(s, e) ? s : e, interior); |
| 1205 | } |
| 1206 | |
| 1207 | RS_Vector leftmostArc(const RS_Arc& arc, bool* interior) { |
| 1208 | const RS_Vector c = arc.getCenter(); |
| 1209 | const double r = arc.getRadius(); |
| 1210 | if (!(r > 0.0)) return lexminEndpoints(arc, interior); |
| 1211 | // Walk CCW from a1 by sweep. The angle pi (relative to WCS) corresponds |
| 1212 | // to parameter u in [0, sweep] with (a1 + u) == pi (mod 2*pi) for the |
| 1213 | // non-reversed arc; the reversed arc negates the step but the point of |
| 1214 | // minimum x is unchanged. |
| 1215 | const double a1 = arc.isReversed() ? arc.getAngle2() : arc.getAngle1(); |
| 1216 | const double sweep = std::abs(arc.getAngleLength()); |
| 1217 | if (sweep <= 0.0) return lexminEndpoints(arc, interior); |
| 1218 | double u = std::fmod(kPi - a1, kTwoPi); |
| 1219 | if (u < 0.0) u += kTwoPi; |
| 1220 | const double tolU = 8.0 * kU * std::max(1.0, std::abs(u)); |
| 1221 | if (u > tolU && u < sweep - tolU) { |
| 1222 | if (interior != nullptr) *interior = true; |
| 1223 | return RS_Vector(c.x - r, c.y); |
| 1224 | } |
| 1225 | return lexminEndpoints(arc, interior); |
| 1226 | } |
| 1227 | |
| 1228 | RS_Vector leftmostCircle(const RS_Circle& circ, bool* interior) { |
| 1229 | const RS_Vector c = circ.getCenter(); |
| 1230 | const double r = circ.getRadius(); |
| 1231 | if (!(r > 0.0)) return asVertex(c, interior); |
| 1232 | if (interior != nullptr) *interior = true; |
| 1233 | return RS_Vector(c.x - r, c.y); |
| 1234 | } |
| 1235 | |
| 1236 | RS_Vector ellipsePointAt(const RS_Ellipse& ell, double t) { |
| 1237 | const RS_Vector c = ell.getCenter(); |
| 1238 | const RS_Vector M = ell.getMajorP(); |
| 1239 | const double ratio = ell.getRatio(); |
| 1240 | const RS_Vector N(-ratio * M.y, ratio * M.x); |
| 1241 | return RS_Vector(c.x + M.x * std::cos(t) + N.x * std::sin(t), |
| 1242 | c.y + M.y * std::cos(t) + N.y * std::sin(t)); |
| 1243 | } |
| 1244 | |
| 1245 | RS_Vector leftmostEllipse(const RS_Ellipse& ell, bool* interior) { |
| 1246 | const RS_Vector M = ell.getMajorP(); |
| 1247 | const double ratio = ell.getRatio(); |
| 1248 | const RS_Vector N(-ratio * M.y, ratio * M.x); |
| 1249 | const double Rx = std::hypot(M.x, N.x); |
| 1250 | if (!(Rx > 0.0)) return lexminEndpoints(ell, interior); |
| 1251 | // x(t) = cx + Mx*cos(t) + Nx*sin(t) = cx + Rx*cos(t - beta), beta = atan2(Nx, Mx). |
| 1252 | // Minimum at t = beta + pi, giving x = cx - Rx. |
| 1253 | const double beta = std::atan2(N.x, M.x); |
| 1254 | double tStar = beta + kPi; |
| 1255 | bool fullCircle = !ell.isArc() || |
| 1256 | std::abs(ell.getAngleLength() - kTwoPi) < RS_TOLERANCE_ANGLE1.0e-8; |
| 1257 | if (fullCircle) { |
| 1258 | if (interior != nullptr) *interior = true; |
| 1259 | return ellipsePointAt(ell, tStar); |
| 1260 | } |
| 1261 | const double a1 = ell.isReversed() ? ell.getAngle2() : ell.getAngle1(); |
| 1262 | const double sweep = std::abs(ell.getAngleLength()); |
| 1263 | if (sweep <= 0.0) return lexminEndpoints(ell, interior); |
| 1264 | double u = std::fmod(tStar - a1, kTwoPi); |
| 1265 | if (u < 0.0) u += kTwoPi; |
| 1266 | const double tolU = 8.0 * kU * std::max(1.0, std::abs(u)); |
| 1267 | if (u > tolU && u < sweep - tolU) { |
| 1268 | if (interior != nullptr) *interior = true; |
| 1269 | return ellipsePointAt(ell, a1 + u); |
| 1270 | } |
| 1271 | return lexminEndpoints(ell, interior); |
| 1272 | } |
| 1273 | |
| 1274 | RS_Vector leftmostBezier(const RS_Vector& p0, const RS_Vector& p1, |
| 1275 | const RS_Vector& p2, bool* interior) { |
| 1276 | const double denom = p0.x - 2.0 * p1.x + p2.x; |
| 1277 | if (std::abs(denom) > 1e-14 && denom > 0.0) { |
| 1278 | const double w = (p0.x - p1.x) / denom; |
| 1279 | if (w > 1e-12 && w < 1.0 - 1e-12) { |
| 1280 | const double om = 1.0 - w; |
| 1281 | const RS_Vector p(om * om * p0.x + 2.0 * w * om * p1.x + w * w * p2.x, |
| 1282 | om * om * p0.y + 2.0 * w * om * p1.y + w * w * p2.y); |
| 1283 | if (interior != nullptr) *interior = true; |
| 1284 | return p; |
| 1285 | } |
| 1286 | } |
| 1287 | if (interior != nullptr) *interior = false; |
| 1288 | return lexLess(p0, p2) ? p0 : p2; |
| 1289 | } |
| 1290 | |
| 1291 | RS_Vector leftmostParabola(const LC_Parabola& para, bool* interior) { |
| 1292 | const auto& cp = para.getData().m_controlPoints; |
| 1293 | if (cp.size() < 3) return lexminEndpoints(para, interior); |
| 1294 | return leftmostBezier(cp[0], cp[1], cp[2], interior); |
| 1295 | } |
| 1296 | |
| 1297 | RS_Vector leftmostSplinePoints(const LC_SplinePoints& spline, bool* interior) { |
| 1298 | bool haveBest = false; |
| 1299 | bool bestInterior = false; |
| 1300 | RS_Vector best(false); |
| 1301 | const size_t count = spline.getSegmentCount(); |
| 1302 | for (size_t i = 0; i < count; ++i) { |
| 1303 | LC_SplinePointsSegment seg; |
| 1304 | if (!spline.tryGetSegment(i, seg)) continue; |
| 1305 | bool segInterior = false; |
| 1306 | RS_Vector cand(false); |
| 1307 | switch (seg.kind) { |
| 1308 | case LC_SplinePointsSegment::Kind::Point: |
| 1309 | if (seg.start.valid) cand = seg.start; |
| 1310 | break; |
| 1311 | case LC_SplinePointsSegment::Kind::Line: |
| 1312 | cand = lexLess(seg.start, seg.end) ? seg.start : seg.end; |
| 1313 | break; |
| 1314 | case LC_SplinePointsSegment::Kind::Quadratic: |
| 1315 | cand = leftmostBezier(seg.start, seg.control, seg.end, &segInterior); |
| 1316 | break; |
| 1317 | } |
| 1318 | if (!cand.valid) continue; |
| 1319 | if (!haveBest || lexLess(cand, best)) { |
| 1320 | best = cand; |
| 1321 | bestInterior = segInterior; |
| 1322 | haveBest = true; |
| 1323 | } |
| 1324 | } |
| 1325 | if (!haveBest) return lexminEndpoints(spline, interior); |
| 1326 | if (interior != nullptr) *interior = bestInterior; |
| 1327 | return best; |
| 1328 | } |
| 1329 | |
| 1330 | RS_Vector leftmostHyperbola(const LC_Hyperbola& hyp, bool* interior) { |
| 1331 | if (hyp.isInfinite()) return lexminEndpoints(hyp, interior); |
| 1332 | const auto& d = hyp.getData(); |
| 1333 | const double a = d.majorP.magnitude(); |
| 1334 | const double b = a * d.ratio; |
| 1335 | if (!(a > 0.0) || !(b > 0.0)) return lexminEndpoints(hyp, interior); |
| 1336 | const double th = d.majorP.angle(); |
| 1337 | const double sx = d.reversed ? -1.0 : 1.0; |
| 1338 | // x(phi) = cx + sx*a*cosh(phi)*cos(th) - b*sinh(phi)*sin(th) |
| 1339 | // = cx + Ax*cosh(phi) + Bx*sinh(phi), Ax = sx*a*cos(th), Bx = -b*sin(th) |
| 1340 | // Minimum iff |Bx| < |Ax| and Ax > 0, at tanh(phi*) = -Bx/Ax. |
| 1341 | const double Ax = sx * a * std::cos(th); |
| 1342 | const double Bx = -b * std::sin(th); |
| 1343 | if (Ax > 0.0 && std::abs(Bx) < std::abs(Ax)) { |
| 1344 | const double r = -Bx / Ax; |
| 1345 | if (r > -1.0 && r < 1.0) { |
| 1346 | const double phiStar = std::atanh(r); |
| 1347 | const double phi1 = std::min(d.angle1, d.angle2); |
| 1348 | const double phi2 = std::max(d.angle1, d.angle2); |
| 1349 | if (phiStar > phi1 + 1e-12 && phiStar < phi2 - 1e-12) { |
| 1350 | const RS_Vector p = hyp.getPoint(phiStar, d.reversed); |
| 1351 | if (p.valid) { |
| 1352 | if (interior != nullptr) *interior = true; |
| 1353 | return p; |
| 1354 | } |
| 1355 | } |
| 1356 | } |
| 1357 | } |
| 1358 | return lexminEndpoints(hyp, interior); |
| 1359 | } |
| 1360 | |
| 1361 | } // namespace |
| 1362 | |
| 1363 | RS_Vector leftmostPoint(const RS_Entity& edge, bool* interior) { |
| 1364 | switch (edge.rtti()) { |
| 1365 | case RS2::EntityLine: |
| 1366 | return lexminEndpoints(edge, interior); |
| 1367 | case RS2::EntityArc: |
| 1368 | return leftmostArc(static_cast<const RS_Arc&>(edge), interior); |
| 1369 | case RS2::EntityCircle: |
| 1370 | return leftmostCircle(static_cast<const RS_Circle&>(edge), interior); |
| 1371 | case RS2::EntityEllipse: |
| 1372 | return leftmostEllipse(static_cast<const RS_Ellipse&>(edge), interior); |
| 1373 | case RS2::EntityParabola: |
| 1374 | return leftmostParabola(static_cast<const LC_Parabola&>(edge), interior); |
| 1375 | case RS2::EntitySplinePoints: |
| 1376 | return leftmostSplinePoints(static_cast<const LC_SplinePoints&>(edge), |
| 1377 | interior); |
| 1378 | case RS2::EntityHyperbola: |
| 1379 | return leftmostHyperbola(static_cast<const LC_Hyperbola&>(edge), |
| 1380 | interior); |
| 1381 | default: |
| 1382 | return lexminEndpoints(edge, interior); |
| 1383 | } |
| 1384 | } |
| 1385 | |
| 1386 | RS_Vector travelDirection(const RS_Entity& edge, bool atStart) { |
| 1387 | auto fromAngle = [](double a) { |
| 1388 | return RS_Vector(std::cos(a), std::sin(a)); |
| 1389 | }; |
| 1390 | switch (edge.rtti()) { |
| 1391 | case RS2::EntityLine: { |
| 1392 | const auto& ln = static_cast<const RS_Line&>(edge); |
| 1393 | const RS_Vector d = ln.getEndpoint() - ln.getStartpoint(); |
| 1394 | const double m = d.magnitude(); |
| 1395 | if (!(m > 0.0)) return RS_Vector(0.0, 0.0); |
| 1396 | return d / m; |
| 1397 | } |
| 1398 | case RS2::EntityArc: |
| 1399 | case RS2::EntityEllipse: { |
| 1400 | const double a = atStart ? edge.getDirection1() |
| 1401 | : edge.getDirection2() + kPi; |
| 1402 | return fromAngle(a); |
| 1403 | } |
| 1404 | case RS2::EntitySplinePoints: { |
| 1405 | const double a = atStart ? edge.getDirection1() |
| 1406 | : edge.getDirection2() + kPi; |
| 1407 | return fromAngle(a); |
| 1408 | } |
| 1409 | case RS2::EntityParabola: { |
| 1410 | const double a = atStart ? edge.getDirection1() |
| 1411 | : edge.getDirection2(); |
| 1412 | return fromAngle(a); |
| 1413 | } |
| 1414 | case RS2::EntityHyperbola: { |
| 1415 | const auto& hyp = static_cast<const LC_Hyperbola&>(edge); |
| 1416 | const auto& d = hyp.getData(); |
| 1417 | const double param = atStart ? d.angle1 : d.angle2; |
| 1418 | RS_Vector t = hyp.getTangentDirectionParam(param); |
| 1419 | if (!t.valid) return RS_Vector(0.0, 0.0); |
| 1420 | const double sign = (d.angle2 >= d.angle1) ? 1.0 : -1.0; |
| 1421 | return t * sign; |
| 1422 | } |
| 1423 | case RS2::EntityCircle: |
| 1424 | // A circle has no well-defined start/end travel direction. |
| 1425 | return RS_Vector(0.0, 0.0); |
| 1426 | default: |
| 1427 | return RS_Vector(0.0, 0.0); |
| 1428 | } |
| 1429 | } |
| 1430 | |
| 1431 | double signedCurvature(const RS_Entity& edge, bool atStart) { |
| 1432 | (void)atStart; // signed curvature is uniform along the pieces we handle. |
| 1433 | switch (edge.rtti()) { |
| 1434 | case RS2::EntityLine: |
| 1435 | return 0.0; |
| 1436 | case RS2::EntityArc: { |
| 1437 | const auto& arc = static_cast<const RS_Arc&>(edge); |
| 1438 | const double r = arc.getRadius(); |
| 1439 | if (!(r > 0.0)) return 0.0; |
| 1440 | return arc.isReversed() ? -1.0 / r : 1.0 / r; |
| 1441 | } |
| 1442 | case RS2::EntityCircle: { |
| 1443 | const auto& circ = static_cast<const RS_Circle&>(edge); |
| 1444 | const double r = circ.getRadius(); |
| 1445 | if (!(r > 0.0)) return 0.0; |
| 1446 | return 1.0 / r; |
| 1447 | } |
| 1448 | case RS2::EntityEllipse: { |
| 1449 | const auto& ell = static_cast<const RS_Ellipse&>(edge); |
| 1450 | const RS_Vector M = ell.getMajorP(); |
| 1451 | const double ratio = ell.getRatio(); |
| 1452 | const double a = M.magnitude(); |
| 1453 | const double b = a * ratio; |
| 1454 | if (!(a > 0.0) || !(b > 0.0)) return 0.0; |
| 1455 | const double k = b / (a * a); |
| 1456 | return ell.isReversed() ? -k : k; |
| 1457 | } |
| 1458 | case RS2::EntitySplinePoints: |
| 1459 | case RS2::EntityParabola: |
| 1460 | case RS2::EntityHyperbola: |
| 1461 | default: |
| 1462 | return 0.0; |
| 1463 | } |
| 1464 | } |
| 1465 | |
| 1466 | } // namespace LC_LoopUtils |