gfx/2d/BezierUtils.cpp

   1 /* -*- Mode: C++; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
   2 /* vim: set ts=8 sts=2 et sw=2 tw=80: */
   3 /* This Source Code Form is subject to the terms of the Mozilla Public
   4  * License, v. 2.0. If a copy of the MPL was not distributed with this
   5  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
   6
   7 #include "BezierUtils.h"
   8
   9 #include "PathHelpers.h"
  10
  11 namespace mozilla {
  12 namespace gfx {
  13
  14 Point GetBezierPoint(const Bezier& aBezier, Float t) {
  15   Float s = 1.0f - t;
  16
  17   return Point(aBezier.mPoints[0].x * s * s * s +
  18                    3.0f * aBezier.mPoints[1].x * t * s * s +
  19                    3.0f * aBezier.mPoints[2].x * t * t * s +
  20                    aBezier.mPoints[3].x * t * t * t,
  21                aBezier.mPoints[0].y * s * s * s +
  22                    3.0f * aBezier.mPoints[1].y * t * s * s +
  23                    3.0f * aBezier.mPoints[2].y * t * t * s +
  24                    aBezier.mPoints[3].y * t * t * t);
  25 }
  26
  27 Point GetBezierDifferential(const Bezier& aBezier, Float t) {
  28   // Return P'(t).
  29
  30   Float s = 1.0f - t;
  31
  32   return Point(
  33       -3.0f * ((aBezier.mPoints[0].x - aBezier.mPoints[1].x) * s * s +
  34                2.0f * (aBezier.mPoints[1].x - aBezier.mPoints[2].x) * t * s +
  35                (aBezier.mPoints[2].x - aBezier.mPoints[3].x) * t * t),
  36       -3.0f * ((aBezier.mPoints[0].y - aBezier.mPoints[1].y) * s * s +
  37                2.0f * (aBezier.mPoints[1].y - aBezier.mPoints[2].y) * t * s +
  38                (aBezier.mPoints[2].y - aBezier.mPoints[3].y) * t * t));
  39 }
  40
  41 Point GetBezierDifferential2(const Bezier& aBezier, Float t) {
  42   // Return P''(t).
  43
  44   Float s = 1.0f - t;
  45
  46   return Point(6.0f * ((aBezier.mPoints[0].x - aBezier.mPoints[1].x) * s -
  47                        (aBezier.mPoints[1].x - aBezier.mPoints[2].x) * (s - t) -
  48                        (aBezier.mPoints[2].x - aBezier.mPoints[3].x) * t),
  49                6.0f * ((aBezier.mPoints[0].y - aBezier.mPoints[1].y) * s -
  50                        (aBezier.mPoints[1].y - aBezier.mPoints[2].y) * (s - t) -
  51                        (aBezier.mPoints[2].y - aBezier.mPoints[3].y) * t));
  52 }
  53
  54 Float GetBezierLength(const Bezier& aBezier, Float a, Float b) {
  55   if (a < 0.5f && b > 0.5f) {
  56     // To increase the accuracy, split into two parts.
  57     return GetBezierLength(aBezier, a, 0.5f) +
  58            GetBezierLength(aBezier, 0.5f, b);
  59   }
  60
  61   // Calculate length of simple bezier curve with Simpson's rule.
  62   //            _
  63   //           /  b
  64   // length =  |    |P'(x)| dx
  65   //          _/  a
  66   //
  67   //          b - a                   a + b
  68   //        = ----- [ |P'(a)| + 4 |P'(-----)| + |P'(b)| ]
  69   //            6                       2
  70
  71   Float fa = GetBezierDifferential(aBezier, a).Length();
  72   Float fab = GetBezierDifferential(aBezier, (a + b) / 2.0f).Length();
  73   Float fb = GetBezierDifferential(aBezier, b).Length();
  74
  75   return (b - a) / 6.0f * (fa + 4.0f * fab + fb);
  76 }
  77
  78 static void SplitBezierA(Bezier* aSubBezier, const Bezier& aBezier, Float t) {
  79   // Split bezier curve into [0,t] and [t,1] parts, and return [0,t] part.
  80
  81   Float s = 1.0f - t;
  82
  83   Point tmp1;
  84   Point tmp2;
  85
  86   aSubBezier->mPoints[0] = aBezier.mPoints[0];
  87
  88   aSubBezier->mPoints[1] = aBezier.mPoints[0] * s + aBezier.mPoints[1] * t;
  89   tmp1 = aBezier.mPoints[1] * s + aBezier.mPoints[2] * t;
  90   tmp2 = aBezier.mPoints[2] * s + aBezier.mPoints[3] * t;
  91
  92   aSubBezier->mPoints[2] = aSubBezier->mPoints[1] * s + tmp1 * t;
  93   tmp1 = tmp1 * s + tmp2 * t;
  94
  95   aSubBezier->mPoints[3] = aSubBezier->mPoints[2] * s + tmp1 * t;
  96 }
  97
  98 static void SplitBezierB(Bezier* aSubBezier, const Bezier& aBezier, Float t) {
  99   // Split bezier curve into [0,t] and [t,1] parts, and return [t,1] part.
 100
 101   Float s = 1.0f - t;
 102
 103   Point tmp1;
 104   Point tmp2;
 105
 106   aSubBezier->mPoints[3] = aBezier.mPoints[3];
 107
 108   aSubBezier->mPoints[2] = aBezier.mPoints[2] * s + aBezier.mPoints[3] * t;
 109   tmp1 = aBezier.mPoints[1] * s + aBezier.mPoints[2] * t;
 110   tmp2 = aBezier.mPoints[0] * s + aBezier.mPoints[1] * t;
 111
 112   aSubBezier->mPoints[1] = tmp1 * s + aSubBezier->mPoints[2] * t;
 113   tmp1 = tmp2 * s + tmp1 * t;
 114
 115   aSubBezier->mPoints[0] = tmp1 * s + aSubBezier->mPoints[1] * t;
 116 }
 117
 118 void GetSubBezier(Bezier* aSubBezier, const Bezier& aBezier, Float t1,
 119                   Float t2) {
 120   Bezier tmp;
 121   SplitBezierB(&tmp, aBezier, t1);
 122
 123   Float range = 1.0f - t1;
 124   if (range == 0.0f) {
 125     *aSubBezier = tmp;
 126   } else {
 127     SplitBezierA(aSubBezier, tmp, (t2 - t1) / range);
 128   }
 129 }
 130
 131 static Point BisectBezierNearestPoint(const Bezier& aBezier,
 132                                       const Point& aTarget, Float* aT) {
 133   // Find a nearest point on bezier curve with Binary search.
 134   // Called from FindBezierNearestPoint.
 135
 136   Float lower = 0.0f;
 137   Float upper = 1.0f;
 138   Float t;
 139
 140   Point P, lastP;
 141   const size_t MAX_LOOP = 32;
 142   const Float DIST_MARGIN = 0.1f;
 143   const Float DIST_MARGIN_SQUARE = DIST_MARGIN * DIST_MARGIN;
 144   const Float DIFF = 0.0001f;
 145   for (size_t i = 0; i < MAX_LOOP; i++) {
 146     t = (upper + lower) / 2.0f;
 147     P = GetBezierPoint(aBezier, t);
 148
 149     // Check if it converged.
 150     if (i > 0 && (lastP - P).LengthSquare() < DIST_MARGIN_SQUARE) {
 151       break;
 152     }
 153
 154     Float distSquare = (P - aTarget).LengthSquare();
 155     if ((GetBezierPoint(aBezier, t + DIFF) - aTarget).LengthSquare() <
 156         distSquare) {
 157       lower = t;
 158     } else if ((GetBezierPoint(aBezier, t - DIFF) - aTarget).LengthSquare() <
 159                distSquare) {
 160       upper = t;
 161     } else {
 162       break;
 163     }
 164
 165     lastP = P;
 166   }
 167
 168   if (aT) {
 169     *aT = t;
 170   }
 171
 172   return P;
 173 }
 174
 175 Point FindBezierNearestPoint(const Bezier& aBezier, const Point& aTarget,
 176                              Float aInitialT, Float* aT) {
 177   // Find a nearest point on bezier curve with Newton's method.
 178   // It converges within 4 iterations in most cases.
 179   //
 180   //                   f(t_n)
 181   //  t_{n+1} = t_n - ---------
 182   //                   f'(t_n)
 183   //
 184   //             d                     2
 185   //     f(t) = ---- | P(t) - aTarget |
 186   //             dt
 187
 188   Float t = aInitialT;
 189   Point P;
 190   Point lastP = GetBezierPoint(aBezier, t);
 191
 192   const size_t MAX_LOOP = 4;
 193   const Float DIST_MARGIN = 0.1f;
 194   const Float DIST_MARGIN_SQUARE = DIST_MARGIN * DIST_MARGIN;
 195   for (size_t i = 0; i <= MAX_LOOP; i++) {
 196     Point dP = GetBezierDifferential(aBezier, t);
 197     Point ddP = GetBezierDifferential2(aBezier, t);
 198     Float f = 2.0f * (lastP.DotProduct(dP) - aTarget.DotProduct(dP));
 199     Float df = 2.0f * (dP.DotProduct(dP) + lastP.DotProduct(ddP) -
 200                        aTarget.DotProduct(ddP));
 201     t = t - f / df;
 202     P = GetBezierPoint(aBezier, t);
 203     if ((P - lastP).LengthSquare() < DIST_MARGIN_SQUARE) {
 204       break;
 205     }
 206     lastP = P;
 207
 208     if (i == MAX_LOOP) {
 209       // If aInitialT is too bad, it won't converge in a few iterations,
 210       // fallback to binary search.
 211       return BisectBezierNearestPoint(aBezier, aTarget, aT);
 212     }
 213   }
 214
 215   if (aT) {
 216     *aT = t;
 217   }
 218
 219   return P;
 220 }
 221
 222 void GetBezierPointsForCorner(Bezier* aBezier, Corner aCorner,
 223                               const Point& aCornerPoint,
 224                               const Size& aCornerSize) {
 225   // Calculate bezier control points for elliptic arc.
 226
 227   const Float signsList[4][2] = {
 228       {+1.0f, +1.0f}, {-1.0f, +1.0f}, {-1.0f, -1.0f}, {+1.0f, -1.0f}};
 229   const Float(&signs)[2] = signsList[aCorner];
 230
 231   aBezier->mPoints[0] = aCornerPoint;
 232   aBezier->mPoints[0].x += signs[0] * aCornerSize.width;
 233
 234   aBezier->mPoints[1] = aBezier->mPoints[0];
 235   aBezier->mPoints[1].x -= signs[0] * aCornerSize.width * kKappaFactor;
 236
 237   aBezier->mPoints[3] = aCornerPoint;
 238   aBezier->mPoints[3].y += signs[1] * aCornerSize.height;
 239
 240   aBezier->mPoints[2] = aBezier->mPoints[3];
 241   aBezier->mPoints[2].y -= signs[1] * aCornerSize.height * kKappaFactor;
 242 }
 243
 244 Float GetQuarterEllipticArcLength(Float a, Float b) {
 245   // Calculate the approximate length of a quarter elliptic arc formed by radii
 246   // (a, b), by Ramanujan's approximation of the perimeter p of an ellipse.
 247   //           _                                                     _
 248   //          |                                      2                |
 249   //          |                           3 * (a - b)                 |
 250   //  p =  PI | (a + b) + ------------------------------------------- |
 251   //          |                                 2                 2   |
 252   //          |_           10 * (a + b) + sqrt(a  + 14 * a * b + b ) _|
 253   //
 254   //           _                                                            _
 255   //          |                                         2                    |
 256   //          |                              3 * (a - b)                     |
 257   //    =  PI | (a + b) + -------------------------------------------------- |
 258   //          |                                           2              2   |
 259   //          |_           10 * (a + b) + sqrt(4 * (a + b)  - 3 * (a - b) ) _|
 260   //
 261   //           _                                          _
 262   //          |                          2                 |
 263   //          |                     3 * S                  |
 264   //    =  PI | A + -------------------------------------- |
 265   //          |                               2        2   |
 266   //          |_           10 * A + sqrt(4 * A  - 3 * S ) _|
 267   //
 268   // where A = a + b, S = a - b
 269
 270   Float A = a + b, S = a - b;
 271   Float A2 = A * A, S2 = S * S;
 272   Float p = M_PI * (A + 3.0f * S2 / (10.0f * A + sqrt(4.0f * A2 - 3.0f * S2)));
 273   return p / 4.0f;
 274 }
 275
 276 Float CalculateDistanceToEllipticArc(const Point& P, const Point& normal,
 277                                      const Point& origin, Float width,
 278                                      Float height) {
 279   // Solve following equations with n and return smaller n.
 280   //
 281   // /  (x, y) = P + n * normal
 282   // |
 283   // <   _            _ 2   _            _ 2
 284   // |  | x - origin.x |   | y - origin.y |
 285   // |  | ------------ | + | ------------ | = 1
 286   // \  |_   width    _|   |_   height   _|
 287
 288   Float a = (P.x - origin.x) / width;
 289   Float b = normal.x / width;
 290   Float c = (P.y - origin.y) / height;
 291   Float d = normal.y / height;
 292
 293   Float A = b * b + d * d;
 294   // In the quadratic formulat B would be 2*(a*b+c*d), however we factor the 2
 295   // out Here which cancels out later.
 296   Float B = a * b + c * d;
 297   Float C = a * a + c * c - 1.0;
 298
 299   Float signB = 1.0;
 300   if (B < 0.0) {
 301     signB = -1.0;
 302   }
 303
 304   // 2nd degree polynomials are typically computed using the formulae
 305   // r1 = -(B - sqrt(delta)) / (2 * A)
 306   // r2 = -(B + sqrt(delta)) / (2 * A)
 307   // However B - sqrt(delta) can be an inportant source of precision loss for
 308   // one of the roots when computing the difference between two similar and
 309   // large numbers. To avoid that we pick the root with no precision loss in r1
 310   // and compute r2 using the Citardauq formula.
 311   // Factoring out 2 from B earlier let
 312   Float S = B + signB * sqrt(B * B - A * C);
 313   Float r1 = -S / A;
 314   Float r2 = -C / S;
 315
 316 #ifdef DEBUG
 317   Float epsilon = (Float)0.001;
 318   MOZ_ASSERT(r1 >= -epsilon);
 319   MOZ_ASSERT(r2 >= -epsilon);
 320 #endif
 321
 322   return std::max((r1 < r2 ? r1 : r2), (Float)0.0);
 323 }
 324
 325 }  // namespace gfx
 326 }  // namespace mozilla