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