BezierUtils.cpp (11045B)
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