From 07393cab57ce74a4aae89a31fae9aaa9780fc19d Mon Sep 17 00:00:00 2001 From: "caryclark@google.com" Date: Mon, 8 Apr 2013 11:47:37 +0000 Subject: Add base types for path ops Paths contain lines, quads, and cubics, which are collectively curves. To work with path intersections, intermediary curves are constructed. For now, those intermediates use doubles to guarantee sufficient precision. The DVector, DPoint, DLine, DQuad, and DCubic structs encapsulate these intermediate curves. The DRect and DTriangle structs are created to describe intersectable areas of interest. The Bounds struct inherits from SkRect to create a SkScalar-based rectangle that intersects shared edges. This also includes common math equalities and debugging that the remainder of path ops builds on, as well as a temporary top-level interface in include/pathops/SkPathOps.h. Review URL: https://codereview.chromium.org/12827020 git-svn-id: http://skia.googlecode.com/svn/trunk@8551 2bbb7eff-a529-9590-31e7-b0007b416f81 --- src/pathops/SkQuarticRoot.cpp | 165 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 src/pathops/SkQuarticRoot.cpp (limited to 'src/pathops/SkQuarticRoot.cpp') diff --git a/src/pathops/SkQuarticRoot.cpp b/src/pathops/SkQuarticRoot.cpp new file mode 100644 index 0000000000..09a92c66a2 --- /dev/null +++ b/src/pathops/SkQuarticRoot.cpp @@ -0,0 +1,165 @@ +// from http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c +/* + * Roots3And4.c + * + * Utility functions to find cubic and quartic roots, + * coefficients are passed like this: + * + * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0 + * + * The functions return the number of non-complex roots and + * put the values into the s array. + * + * Author: Jochen Schwarze (schwarze@isa.de) + * + * Jan 26, 1990 Version for Graphics Gems + * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic + * (reported by Mark Podlipec), + * Old-style function definitions, + * IsZero() as a macro + * Nov 23, 1990 Some systems do not declare acos() and cbrt() in + * , though the functions exist in the library. + * If large coefficients are used, EQN_EPS should be + * reduced considerably (e.g. to 1E-30), results will be + * correct but multiple roots might be reported more + * than once. + */ + +#include "SkPathOpsCubic.h" +#include "SkPathOpsQuad.h" +#include "SkQuarticRoot.h" + +int SkReducedQuarticRoots(const double t4, const double t3, const double t2, const double t1, + const double t0, const bool oneHint, double roots[4]) { +#ifdef SK_DEBUG + // create a string mathematica understands + // GDB set print repe 15 # if repeated digits is a bother + // set print elements 400 # if line doesn't fit + char str[1024]; + sk_bzero(str, sizeof(str)); + SK_SNPRINTF(str, sizeof(str), + "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", + t4, t3, t2, t1, t0); + mathematica_ize(str, sizeof(str)); +#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA + SkDebugf("%s\n", str); +#endif +#endif + if (approximately_zero_when_compared_to(t4, t0) // 0 is one root + && approximately_zero_when_compared_to(t4, t1) + && approximately_zero_when_compared_to(t4, t2)) { + if (approximately_zero_when_compared_to(t3, t0) + && approximately_zero_when_compared_to(t3, t1) + && approximately_zero_when_compared_to(t3, t2)) { + return SkDQuad::RootsReal(t2, t1, t0, roots); + } + if (approximately_zero_when_compared_to(t4, t3)) { + return SkDCubic::RootsReal(t3, t2, t1, t0, roots); + } + } + if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1)) // 0 is one root + // && approximately_zero_when_compared_to(t0, t2) + && approximately_zero_when_compared_to(t0, t3) + && approximately_zero_when_compared_to(t0, t4)) { + int num = SkDCubic::RootsReal(t4, t3, t2, t1, roots); + for (int i = 0; i < num; ++i) { + if (approximately_zero(roots[i])) { + return num; + } + } + roots[num++] = 0; + return num; + } + if (oneHint) { + SkASSERT(approximately_zero(t4 + t3 + t2 + t1 + t0)); // 1 is one root + // note that -C == A + B + D + E + int num = SkDCubic::RootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); + for (int i = 0; i < num; ++i) { + if (approximately_equal(roots[i], 1)) { + return num; + } + } + roots[num++] = 1; + return num; + } + return -1; +} + +int SkQuarticRootsReal(int firstCubicRoot, const double A, const double B, const double C, + const double D, const double E, double s[4]) { + double u, v; + /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ + const double invA = 1 / A; + const double a = B * invA; + const double b = C * invA; + const double c = D * invA; + const double d = E * invA; + /* substitute x = y - a/4 to eliminate cubic term: + x^4 + px^2 + qx + r = 0 */ + const double a2 = a * a; + const double p = -3 * a2 / 8 + b; + const double q = a2 * a / 8 - a * b / 2 + c; + const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d; + int num; + if (approximately_zero(r)) { + /* no absolute term: y(y^3 + py + q) = 0 */ + num = SkDCubic::RootsReal(1, 0, p, q, s); + s[num++] = 0; + } else { + /* solve the resolvent cubic ... */ + double cubicRoots[3]; + int roots = SkDCubic::RootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots); + int index; + /* ... and take one real solution ... */ + double z; + num = 0; + int num2 = 0; + for (index = firstCubicRoot; index < roots; ++index) { + z = cubicRoots[index]; + /* ... to build two quadric equations */ + u = z * z - r; + v = 2 * z - p; + if (approximately_zero_squared(u)) { + u = 0; + } else if (u > 0) { + u = sqrt(u); + } else { + continue; + } + if (approximately_zero_squared(v)) { + v = 0; + } else if (v > 0) { + v = sqrt(v); + } else { + continue; + } + num = SkDQuad::RootsReal(1, q < 0 ? -v : v, z - u, s); + num2 = SkDQuad::RootsReal(1, q < 0 ? v : -v, z + u, s + num); + if (!((num | num2) & 1)) { + break; // prefer solutions without single quad roots + } + } + num += num2; + if (!num) { + return 0; // no valid cubic root + } + } + /* resubstitute */ + const double sub = a / 4; + for (int i = 0; i < num; ++i) { + s[i] -= sub; + } + // eliminate duplicates + for (int i = 0; i < num - 1; ++i) { + for (int j = i + 1; j < num; ) { + if (AlmostEqualUlps(s[i], s[j])) { + if (j < --num) { + s[j] = s[num]; + } + } else { + ++j; + } + } + } + return num; +} -- cgit v1.2.3