1/* 2 * Copyright 2016 Google Inc. 3 * 4 * Use of this source code is governed by a BSD-style license that can be 5 * found in the LICENSE file. 6 */ 7 8#include "SkCurveMeasure.h" 9#include "SkGeometry.h" 10 11// for abs 12#include <cmath> 13 14#define UNIMPLEMENTED SkDEBUGF(("%s:%d unimplemented\n", __FILE__, __LINE__)) 15 16/// Used inside SkCurveMeasure::getTime's Newton's iteration 17static inline SkPoint evaluate(const SkPoint pts[4], SkSegType segType, 18 SkScalar t) { 19 SkPoint pos; 20 switch (segType) { 21 case kQuad_SegType: 22 pos = SkEvalQuadAt(pts, t); 23 break; 24 case kLine_SegType: 25 pos = SkPoint::Make(SkScalarInterp(pts[0].x(), pts[1].x(), t), 26 SkScalarInterp(pts[0].y(), pts[1].y(), t)); 27 break; 28 case kCubic_SegType: 29 SkEvalCubicAt(pts, t, &pos, nullptr, nullptr); 30 break; 31 case kConic_SegType: { 32 SkConic conic(pts, pts[3].x()); 33 conic.evalAt(t, &pos); 34 } 35 break; 36 default: 37 UNIMPLEMENTED; 38 } 39 40 return pos; 41} 42 43/// Used inside SkCurveMeasure::getTime's Newton's iteration 44static inline SkVector evaluateDerivative(const SkPoint pts[4], 45 SkSegType segType, SkScalar t) { 46 SkVector tan; 47 switch (segType) { 48 case kQuad_SegType: 49 tan = SkEvalQuadTangentAt(pts, t); 50 break; 51 case kLine_SegType: 52 tan = pts[1] - pts[0]; 53 break; 54 case kCubic_SegType: 55 SkEvalCubicAt(pts, t, nullptr, &tan, nullptr); 56 break; 57 case kConic_SegType: { 58 SkConic conic(pts, pts[3].x()); 59 conic.evalAt(t, nullptr, &tan); 60 } 61 break; 62 default: 63 UNIMPLEMENTED; 64 } 65 66 return tan; 67} 68/// Used in ArcLengthIntegrator::computeLength 69static inline Sk8f evaluateDerivativeLength(const Sk8f& ts, 70 const float (&xCoeff)[3][8], 71 const float (&yCoeff)[3][8], 72 const SkSegType segType) { 73 Sk8f x; 74 Sk8f y; 75 76 Sk8f x0 = Sk8f::Load(&xCoeff[0]), 77 x1 = Sk8f::Load(&xCoeff[1]), 78 x2 = Sk8f::Load(&xCoeff[2]); 79 80 Sk8f y0 = Sk8f::Load(&yCoeff[0]), 81 y1 = Sk8f::Load(&yCoeff[1]), 82 y2 = Sk8f::Load(&yCoeff[2]); 83 84 switch (segType) { 85 case kQuad_SegType: 86 x = x0*ts + x1; 87 y = y0*ts + y1; 88 break; 89 case kCubic_SegType: 90 x = (x0*ts + x1)*ts + x2; 91 y = (y0*ts + y1)*ts + y2; 92 break; 93 case kConic_SegType: 94 UNIMPLEMENTED; 95 break; 96 default: 97 UNIMPLEMENTED; 98 } 99 100 x = x * x; 101 y = y * y; 102 103 return (x + y).sqrt(); 104} 105 106ArcLengthIntegrator::ArcLengthIntegrator(const SkPoint* pts, SkSegType segType) 107 : fSegType(segType) { 108 switch (fSegType) { 109 case kQuad_SegType: { 110 float Ax = pts[0].x(); 111 float Bx = pts[1].x(); 112 float Cx = pts[2].x(); 113 float Ay = pts[0].y(); 114 float By = pts[1].y(); 115 float Cy = pts[2].y(); 116 117 // precompute coefficients for derivative 118 Sk8f(2*(Ax - 2*Bx + Cx)).store(&xCoeff[0]); 119 Sk8f(2*(Bx - Ax)).store(&xCoeff[1]); 120 121 Sk8f(2*(Ay - 2*By + Cy)).store(&yCoeff[0]); 122 Sk8f(2*(By - Ay)).store(&yCoeff[1]); 123 } 124 break; 125 case kCubic_SegType: 126 { 127 float Ax = pts[0].x(); 128 float Bx = pts[1].x(); 129 float Cx = pts[2].x(); 130 float Dx = pts[3].x(); 131 float Ay = pts[0].y(); 132 float By = pts[1].y(); 133 float Cy = pts[2].y(); 134 float Dy = pts[3].y(); 135 136 // precompute coefficients for derivative 137 Sk8f(3*(-Ax + 3*(Bx - Cx) + Dx)).store(&xCoeff[0]); 138 Sk8f(6*(Ax - 2*Bx + Cx)).store(&xCoeff[1]); 139 Sk8f(3*(-Ax + Bx)).store(&xCoeff[2]); 140 141 Sk8f(3*(-Ay + 3*(By - Cy) + Dy)).store(&yCoeff[0]); 142 Sk8f(6*(Ay - 2*By + Cy)).store(&yCoeff[1]); 143 Sk8f(3*(-Ay + By)).store(&yCoeff[2]); 144 } 145 break; 146 case kConic_SegType: 147 UNIMPLEMENTED; 148 break; 149 default: 150 UNIMPLEMENTED; 151 } 152} 153 154// We use Gaussian quadrature 155// (https://en.wikipedia.org/wiki/Gaussian_quadrature) 156// to approximate the arc length integral here, because it is amenable to SIMD. 157SkScalar ArcLengthIntegrator::computeLength(SkScalar t) { 158 SkScalar length = 0.0f; 159 160 Sk8f lengths = evaluateDerivativeLength(absc*t, xCoeff, yCoeff, fSegType); 161 lengths = weights*lengths; 162 // is it faster or more accurate to sum and then multiply or vice versa? 163 lengths = lengths*(t*0.5f); 164 165 // Why does SkNx index with ints? does negative index mean something? 166 for (int i = 0; i < 8; i++) { 167 length += lengths[i]; 168 } 169 return length; 170} 171 172SkCurveMeasure::SkCurveMeasure(const SkPoint* pts, SkSegType segType) 173 : fSegType(segType) { 174 switch (fSegType) { 175 case SkSegType::kQuad_SegType: 176 for (size_t i = 0; i < 3; i++) { 177 fPts[i] = pts[i]; 178 } 179 break; 180 case SkSegType::kLine_SegType: 181 fPts[0] = pts[0]; 182 fPts[1] = pts[1]; 183 fLength = (fPts[1] - fPts[0]).length(); 184 break; 185 case SkSegType::kCubic_SegType: 186 for (size_t i = 0; i < 4; i++) { 187 fPts[i] = pts[i]; 188 } 189 break; 190 case SkSegType::kConic_SegType: 191 for (size_t i = 0; i < 4; i++) { 192 fPts[i] = pts[i]; 193 } 194 break; 195 default: 196 UNIMPLEMENTED; 197 break; 198 } 199 if (kLine_SegType != segType) { 200 fIntegrator = ArcLengthIntegrator(fPts, fSegType); 201 } 202} 203 204SkScalar SkCurveMeasure::getLength() { 205 if (-1.0f == fLength) { 206 fLength = fIntegrator.computeLength(1.0f); 207 } 208 return fLength; 209} 210 211// Given an arc length targetLength, we want to determine what t 212// gives us the corresponding arc length along the curve. 213// We do this by letting the arc length integral := f(t) and 214// solving for the root of the equation f(t) - targetLength = 0 215// using Newton's method and lerp-bisection. 216// The computationally expensive parts are the integral approximation 217// at each step, and computing the derivative of the arc length integral, 218// which is equal to the length of the tangent (so we have to do a sqrt). 219 220SkScalar SkCurveMeasure::getTime(SkScalar targetLength) { 221 if (targetLength <= 0.0f) { 222 return 0.0f; 223 } 224 225 SkScalar currentLength = getLength(); 226 227 if (targetLength > currentLength || (SkScalarNearlyEqual(targetLength, currentLength))) { 228 return 1.0f; 229 } 230 if (kLine_SegType == fSegType) { 231 return targetLength / currentLength; 232 } 233 234 // initial estimate of t is percentage of total length 235 SkScalar currentT = targetLength / currentLength; 236 SkScalar prevT = -1.0f; 237 SkScalar newT; 238 239 SkScalar minT = 0.0f; 240 SkScalar maxT = 1.0f; 241 242 int iterations = 0; 243 while (iterations < kNewtonIters + kBisectIters) { 244 currentLength = fIntegrator.computeLength(currentT); 245 SkScalar lengthDiff = currentLength - targetLength; 246 247 // Update root bounds. 248 // If lengthDiff is positive, we have overshot the target, so 249 // we know the current t is an upper bound, and similarly 250 // for the lower bound. 251 if (lengthDiff > 0.0f) { 252 if (currentT < maxT) { 253 maxT = currentT; 254 } 255 } else { 256 if (currentT > minT) { 257 minT = currentT; 258 } 259 } 260 261 // We have a tolerance on both the absolute value of the difference and 262 // on the t value 263 // because we may not have enough precision in the t to get close enough 264 // in the length. 265 if ((std::abs(lengthDiff) < kTolerance) || 266 (std::abs(prevT - currentT) < kTolerance)) { 267 break; 268 } 269 270 prevT = currentT; 271 if (iterations < kNewtonIters) { 272 // This is just newton's formula. 273 SkScalar dt = evaluateDerivative(fPts, fSegType, currentT).length(); 274 newT = currentT - (lengthDiff / dt); 275 276 // If newT is out of bounds, bisect inside newton. 277 if ((newT < 0.0f) || (newT > 1.0f)) { 278 newT = (minT + maxT) * 0.5f; 279 } 280 } else if (iterations < kNewtonIters + kBisectIters) { 281 if (lengthDiff > 0.0f) { 282 maxT = currentT; 283 } else { 284 minT = currentT; 285 } 286 // TODO(hstern) do a lerp here instead of a bisection 287 newT = (minT + maxT) * 0.5f; 288 } else { 289 SkDEBUGF(("%.7f %.7f didn't get close enough after bisection.\n", 290 currentT, currentLength)); 291 break; 292 } 293 currentT = newT; 294 295 SkASSERT(minT <= maxT); 296 297 iterations++; 298 } 299 300 // debug. is there an SKDEBUG or something for ifdefs? 301 fIters = iterations; 302 303 return currentT; 304} 305 306void SkCurveMeasure::getPosTanTime(SkScalar targetLength, SkPoint* pos, 307 SkVector* tan, SkScalar* time) { 308 SkScalar t = getTime(targetLength); 309 310 if (time) { 311 *time = t; 312 } 313 if (pos) { 314 *pos = evaluate(fPts, fSegType, t); 315 } 316 if (tan) { 317 *tan = evaluateDerivative(fPts, fSegType, t); 318 } 319} 320