src/pathops/SkPathOpsConic.cpp

/*
 * Copyright 2015 Google Inc.
 *
 * Use of this source code is governed by a BSD-style license that can be
 * found in the LICENSE file.
 */
#include "SkIntersections.h"
#include "SkLineParameters.h"
#include "SkPathOpsConic.h"
#include "SkPathOpsCubic.h"
#include "SkPathOpsQuad.h"

// cribbed from the float version in SkGeometry.cpp
static void conic_deriv_coeff(const double src[],
                              SkScalar w,
                              double coeff[3]) {
    const double P20 = src[4] - src[0];
    const double P10 = src[2] - src[0];
    const double wP10 = w * P10;
    coeff[0] = w * P20 - P20;
    coeff[1] = P20 - 2 * wP10;
    coeff[2] = wP10;
}

static double conic_eval_tan(const double coord[], SkScalar w, double t) {
    double coeff[3];
    conic_deriv_coeff(coord, w, coeff);
    return t * (t * coeff[0] + coeff[1]) + coeff[2];
}

int SkDConic::FindExtrema(const double src[], SkScalar w, double t[1]) {
    double coeff[3];
    conic_deriv_coeff(src, w, coeff);

    double tValues[2];
    int roots = SkDQuad::RootsValidT(coeff[0], coeff[1], coeff[2], tValues);
    // In extreme cases, the number of roots returned can be 2. Pathops
    // will fail later on, so there's no advantage to plumbing in an error
    // return here.
    // SkASSERT(0 == roots || 1 == roots);

    if (1 == roots) {
        t[0] = tValues[0];
        return 1;
    }
    return 0;
}

SkDVector SkDConic::dxdyAtT(double t) const {
    SkDVector result = {
        conic_eval_tan(&fPts[0].fX, fWeight, t),
        conic_eval_tan(&fPts[0].fY, fWeight, t)
    };
    if (result.fX == 0 && result.fY == 0) {
        if (zero_or_one(t)) {
            result = fPts[2] - fPts[0];
        } else {
            // incomplete
            SkDebugf("!k");
        }
    }
    return result;
}

static double conic_eval_numerator(const double src[], SkScalar w, double t) {
    SkASSERT(src);
    SkASSERT(t >= 0 && t <= 1);
    double src2w = src[2] * w;
    double C = src[0];
    double A = src[4] - 2 * src2w + C;
    double B = 2 * (src2w - C);
    return (A * t + B) * t + C;
}


static double conic_eval_denominator(SkScalar w, double t) {
    double B = 2 * (w - 1);
    double C = 1;
    double A = -B;
    return (A * t + B) * t + C;
}

bool SkDConic::hullIntersects(const SkDCubic& cubic, bool* isLinear) const {
    return cubic.hullIntersects(*this, isLinear);
}

SkDPoint SkDConic::ptAtT(double t) const {
    if (t == 0) {
        return fPts[0];
    }
    if (t == 1) {
        return fPts[2];
    }
    double denominator = conic_eval_denominator(fWeight, t);
    SkDPoint result = {
        conic_eval_numerator(&fPts[0].fX, fWeight, t) / denominator,
        conic_eval_numerator(&fPts[0].fY, fWeight, t) / denominator
    };
    return result;
}

/* see quad subdivide for point rationale */
/* w rationale : the mid point between t1 and t2 could be determined from the computed a/b/c
   values if the computed w was known. Since we know the mid point at (t1+t2)/2, we'll assume
   that it is the same as the point on the new curve t==(0+1)/2.

    d / dz == conic_poly(dst, unknownW, .5) / conic_weight(unknownW, .5);

    conic_poly(dst, unknownW, .5)
                  =   a / 4 + (b * unknownW) / 2 + c / 4
                  =  (a + c) / 4 + (bx * unknownW) / 2

    conic_weight(unknownW, .5)
                  =   unknownW / 2 + 1 / 2

    d / dz                  == ((a + c) / 2 + b * unknownW) / (unknownW + 1)
    d / dz * (unknownW + 1) ==  (a + c) / 2 + b * unknownW
              unknownW       = ((a + c) / 2 - d / dz) / (d / dz - b)

    Thus, w is the ratio of the distance from the mid of end points to the on-curve point, and the
    distance of the on-curve point to the control point.
 */
SkDConic SkDConic::subDivide(double t1, double t2) const {
    double ax, ay, az;
    if (t1 == 0) {
        ax = fPts[0].fX;
        ay = fPts[0].fY;
        az = 1;
    } else if (t1 != 1) {
        ax = conic_eval_numerator(&fPts[0].fX, fWeight, t1);
        ay = conic_eval_numerator(&fPts[0].fY, fWeight, t1);
        az = conic_eval_denominator(fWeight, t1);
    } else {
        ax = fPts[2].fX;
        ay = fPts[2].fY;
        az = 1;
    }
    double midT = (t1 + t2) / 2;
    double dx = conic_eval_numerator(&fPts[0].fX, fWeight, midT);
    double dy = conic_eval_numerator(&fPts[0].fY, fWeight, midT);
    double dz = conic_eval_denominator(fWeight, midT);
    double cx, cy, cz;
    if (t2 == 1) {
        cx = fPts[2].fX;
        cy = fPts[2].fY;
        cz = 1;
    } else if (t2 != 0) {
        cx = conic_eval_numerator(&fPts[0].fX, fWeight, t2);
        cy = conic_eval_numerator(&fPts[0].fY, fWeight, t2);
        cz = conic_eval_denominator(fWeight, t2);
    } else {
        cx = fPts[0].fX;
        cy = fPts[0].fY;
        cz = 1;
    }
    double bx = 2 * dx - (ax + cx) / 2;
    double by = 2 * dy - (ay + cy) / 2;
    double bz = 2 * dz - (az + cz) / 2;
    SkDConic dst = {{{{ax / az, ay / az}, {bx / bz, by / bz}, {cx / cz, cy / cz}}
            SkDEBUGPARAMS(fPts.fDebugGlobalState) },
            SkDoubleToScalar(bz / sqrt(az * cz)) };
    return dst;
}

SkDPoint SkDConic::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2,
        SkScalar* weight) const {
    SkDConic chopped = this->subDivide(t1, t2);
    *weight = chopped.fWeight;
    return chopped[1];
}