Home | History | Annotate | Download | only in utils 
      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
     17 static 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
     44 static 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
     69 static 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 
    106 ArcLengthIntegrator::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.
    157 SkScalar 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 
    172 SkCurveMeasure::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 
    204 SkScalar 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 
    220 SkScalar 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 
    306 void 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