X-Git-Url: https://git.ucc.asn.au/?p=ipdf%2Fcode.git;a=blobdiff_plain;f=src%2Fbezier.cpp;h=83db943692e39a2f0509425234f2c43e82c49cdd;hp=1c836602ea2e60bbe7e9da455c4e9d282aae0c9f;hb=HEAD;hpb=3172dd5af487e0f8a6e5cd5439dea594b9cbd7c9 diff --git a/src/bezier.cpp b/src/bezier.cpp index 1c83660..83db943 100644 --- a/src/bezier.cpp +++ b/src/bezier.cpp @@ -4,11 +4,108 @@ #include #include + + using namespace std; namespace IPDF { +vector SolveQuadratic(const BReal & a, const BReal & b, const BReal & c, const BReal & min, const BReal & max) +{ + vector roots; roots.reserve(2); + if (a == BReal(0) && b != BReal(0)) + { + roots.push_back(-c/b); + return roots; + } + BReal disc(b*b - BReal(4)*a*c); + if (disc < BReal(0)) + { + return roots; + } + else if (disc == BReal(0)) + { + BReal x(-b/BReal(2)*a); + if (x >= min && x <= max) + roots.push_back(x); + return roots; + } + + BReal x0((-b - Sqrt(b*b - BReal(4)*a*c))/(BReal(2)*a)); + BReal x1((-b + Sqrt(b*b - BReal(4)*a*c))/(BReal(2)*a)); + if (x0 > x1) + { + BReal tmp(x0); + x0 = x1; + x1 = tmp; + } + if (x0 >= min && x0 <= max) + roots.push_back(x0); + if (x1 >= min && x1 <= max) + roots.push_back(x1); + return roots; +} + +/** + * Finds the root (if it exists) in a monotonicly in(de)creasing segment of a Cubic + */ + +static void CubicSolveSegment(vector & roots, const BReal & a, const BReal & b, const BReal & c, const BReal & d, BReal & tl, BReal & tu, const BReal & delta) +{ + BReal l = a*tl*tl*tl + b*tl*tl + c*tl + d; + BReal u = a*tu*tu*tu + b*tu*tu + c*tu + d; + if ((l < BReal(0) && u < BReal(0)) || (l > BReal(0) && u > BReal(0))) + { + //Debug("Discarding segment (no roots) l = %f (%f), u = %f (%f)", Double(tl), Double(l), Double(tu), Double(u)); + //return; + } + + bool negative = (u < l); // lower point > 0, upper point < 0 + //Debug("%ft^3 + %ft^2 + %ft + %f is negative (%f < %f) %d", Double(a),Double(b),Double(c),Double(d),Double(u),Double(l), negative); + while (tu - tl > delta) + { + BReal t(tu+tl); + t /= 2; + BReal m = a*t*t*t + b*t*t + c*t + d; + if (m > BReal(0)) + { + if (negative) + tl = t; + else + tu = t; + } + else if (negative) + { + tu = t; + } + else + { + tl = t; + } + //Debug("Delta is %f (%f - %f -> %f)", tu-tl, tu, tl, t); + } + roots.push_back(tl); +} +vector SolveCubic(const BReal & a, const BReal & b, const BReal & c, const BReal & d, const BReal & min, const BReal & max, const BReal & delta) +{ + vector roots; roots.reserve(3); + BReal tu(max); + BReal tl(min); + vector turns(SolveQuadratic(a*BReal(3), b*BReal(2), c)); + //Debug("%u turning points", turns.size()); + for (unsigned i = 1; i < turns.size(); ++i) + { + //if (tl > max) break; + tu = turns[i]; + CubicSolveSegment(roots, a, b, c, d, tl, tu,delta); + tl = turns[i]; + } + tu = max; + CubicSolveSegment(roots, a, b, c, d, tl, tu,delta); + return roots; +} + /** * Factorial * Use dynamic programming / recursion @@ -41,9 +138,9 @@ int BinomialCoeff(int n, int k) /** * Bernstein Basis Polynomial */ -Real Bernstein(int k, int n, const Real & u) +BReal Bernstein(int k, int n, const BReal & u) { - return Real(BinomialCoeff(n, k)) * Power(u, k) * Power(Real(1.0) - u, n-k); + return BReal(BinomialCoeff(n, k)) * Power(u, k) * Power(BReal(1.0) - u, n-k); } @@ -52,40 +149,42 @@ Real Bernstein(int k, int n, const Real & u) * In one coordinate direction */ -pair BezierTurningPoints(const Real & p0, const Real & p1, const Real & p2, const Real & p3) +pair BezierTurningPoints(const BReal & p0, const BReal & p1, const BReal & p2, const BReal & p3) { // straight line if (p1 == p2 && p2 == p3) { - return pair(0, 1); + return pair(0, 1); } - Real a = (3*(p1-p2) + p3 - p0); - Real b = 2*(p2 - 2*p1 + p0); - Real c = (p1-p0); - if (a == 0) + BReal a = ((p1-p2)*BReal(3) + p3 - p0); + BReal b = (p2 - p1*BReal(2) + p0)*BReal(2); + BReal c = (p1-p0); + if (a == BReal(0)) { - if (b == 0) - return pair(0,1); - Real t = -c/b; - if (t > 1) t = 1; - if (t < 0) t = 0; - return pair(t, t); + if (b == BReal(0)) + return pair(0,1); + BReal t = -c/b; + if (t > BReal(1)) t = 1; + if (t < BReal(0)) t = 0; + return pair(t, t); } //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c)); - if (b*b - 4*a*c < 0) + if (b*b - a*c*BReal(4) < BReal(0)) { //Debug("No real roots"); - return pair(0,1); + return pair(0,1); } - pair tsols = SolveQuadratic(a, b, c); - if (tsols.first > 1) tsols.first = 1; - if (tsols.first < 0) tsols.first = 0; - if (tsols.second > 1) tsols.second = 1; - if (tsols.second < 0) tsols.second = 0; - return tsols; + vector tsols = SolveQuadratic(a, b, c); + if (tsols.size() == 1) + return pair(tsols[0], tsols[0]); + else if (tsols.size() == 0) + return pair(0,1); + + return pair(tsols[0], tsols[1]); + } -inline bool CompRealByPtr(const Real * a, const Real * b) +inline bool CompBRealByPtr(const BReal * a, const BReal * b) { return (*a) < (*b); } @@ -93,20 +192,20 @@ inline bool CompRealByPtr(const Real * a, const Real * b) /** * Get top most *point* on Bezier curve */ -pair Bezier::GetTop() const +pair Bezier::GetTop() const { - pair tsols = BezierTurningPoints(y0,y1,y2,y3); - Real tx0; Real ty0; - Real tx1; Real ty1; + pair tsols = BezierTurningPoints(y0,y1,y2,y3); + BReal tx0; BReal ty0; + BReal tx1; BReal ty1; Evaluate(tx0, ty0, tsols.first); Evaluate(tx1, ty1, tsols.second); - vector v(4); + vector v(4); v[0] = &y0; v[1] = &y3; v[2] = &ty0; v[3] = &ty1; - sort(v.begin(), v.end(), CompRealByPtr); - pair result; + sort(v.begin(), v.end(), CompBRealByPtr); + pair result; result.second = *v[0]; if (v[0] == &y0) { @@ -130,20 +229,20 @@ pair Bezier::GetTop() const /** * Get bottom most *point* on Bezier curve */ -pair Bezier::GetBottom() const +pair Bezier::GetBottom() const { - pair tsols = BezierTurningPoints(y0,y1,y2,y3); - Real tx0; Real ty0; - Real tx1; Real ty1; + pair tsols = BezierTurningPoints(y0,y1,y2,y3); + BReal tx0; BReal ty0; + BReal tx1; BReal ty1; Evaluate(tx0, ty0, tsols.first); Evaluate(tx1, ty1, tsols.second); - vector v(4); + vector v(4); v[0] = &y0; v[1] = &y3; v[2] = &ty0; v[3] = &ty1; - sort(v.begin(), v.end(), CompRealByPtr); - pair result; + sort(v.begin(), v.end(), CompBRealByPtr); + pair result; result.second = *v[3]; if (v[3] == &y0) { @@ -167,20 +266,20 @@ pair Bezier::GetBottom() const /** * Get left most *point* on Bezier curve */ -pair Bezier::GetLeft() const +pair Bezier::GetLeft() const { - pair tsols = BezierTurningPoints(x0,x1,x2,x3); - Real tx0; Real ty0; - Real tx1; Real ty1; + pair tsols = BezierTurningPoints(x0,x1,x2,x3); + BReal tx0; BReal ty0; + BReal tx1; BReal ty1; Evaluate(tx0, ty0, tsols.first); Evaluate(tx1, ty1, tsols.second); - vector v(4); + vector v(4); v[0] = &x0; v[1] = &x3; v[2] = &tx0; v[3] = &tx1; - sort(v.begin(), v.end(), CompRealByPtr); - pair result; + sort(v.begin(), v.end(), CompBRealByPtr); + pair result; result.first = *v[0]; if (v[0] == &x0) { @@ -205,20 +304,20 @@ pair Bezier::GetLeft() const /** * Get left most *point* on Bezier curve */ -pair Bezier::GetRight() const +pair Bezier::GetRight() const { - pair tsols = BezierTurningPoints(x0,x1,x2,x3); - Real tx0; Real ty0; - Real tx1; Real ty1; + pair tsols = BezierTurningPoints(x0,x1,x2,x3); + BReal tx0; BReal ty0; + BReal tx1; BReal ty1; Evaluate(tx0, ty0, tsols.first); Evaluate(tx1, ty1, tsols.second); - vector v(4); + vector v(4); v[0] = &x0; v[1] = &x3; v[2] = &tx0; v[3] = &tx1; - sort(v.begin(), v.end(), CompRealByPtr); - pair result; + sort(v.begin(), v.end(), CompBRealByPtr); + pair result; result.first = *v[3]; if (v[3] == &x0) { @@ -239,28 +338,69 @@ pair Bezier::GetRight() const return result; } +vector Bezier::SolveXParam(const BReal & x) const +{ + BReal d(x0 - x); + BReal c((x1 - x0)*BReal(3)); + BReal b((x2 - x1)*BReal(3) - c); + BReal a(x3 -x0 - c - b); + vector results(SolveCubic(a, b, c, d)); + for (unsigned i = 0; i < results.size(); ++i) + { + Vec2 p; + Evaluate(p.x, p.y, results[i]); + } + return results; +} + + +vector Bezier::SolveYParam(const BReal & y) const +{ + BReal d(y0 - y); + BReal c((y1 - y0)*BReal(3)); + BReal b((y2 - y1)*BReal(3) - c); + BReal a(y3 -y0 - c - b); + vector results(SolveCubic(a, b, c, d)); + for (unsigned i = 0; i < results.size(); ++i) + { + Vec2 p; + Evaluate(p.x, p.y, results[i]); + } + return results; +} + +vector Bezier::Evaluate(const vector & u) const +{ + vector result(u.size()); + for (unsigned i = 0; i < u.size(); ++i) + { + Evaluate(result[i].x, result[i].y, u[i]); + } + return result; +} + /** * Get Bounds Rectangle of Bezier */ -Rect Bezier::SolveBounds() const +BRect Bezier::SolveBounds() const { - Rect result; - pair tsols = BezierTurningPoints(x0, x1, x2, x3); + BRect result; + pair tsols = BezierTurningPoints(x0, x1, x2, x3); - Real tp0; Real tp1; Real o; + BReal tp0; BReal tp1; BReal o; Evaluate(tp0, o, tsols.first); Evaluate(tp1, o, tsols.second); //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1)); - vector v(4); + vector v(4); v[0] = &x0; v[1] = &x3; v[2] = &tp0; v[3] = &tp1; // Not using a lambda to keep this compiling on cabellera - sort(v.begin(), v.end(), CompRealByPtr); + sort(v.begin(), v.end(), CompBRealByPtr); result.x = *(v[0]); result.w = *(v[3]) - result.x; @@ -277,7 +417,7 @@ Rect Bezier::SolveBounds() const v[1] = &y3; v[2] = &tp0; v[3] = &tp1; - sort(v.begin(), v.end(), CompRealByPtr); + sort(v.begin(), v.end(), CompBRealByPtr); result.y = *(v[0]); result.h = *(v[3]) - result.y; @@ -287,3 +427,4 @@ Rect Bezier::SolveBounds() const } } // end namespace +