1/*
2 * Copyright 2015 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#include "SkIntersections.h"
8#include "SkLineParameters.h"
9#include "SkPathOpsConic.h"
10#include "SkPathOpsCubic.h"
11#include "SkPathOpsQuad.h"
12
13// cribbed from the float version in SkGeometry.cpp
14static void conic_deriv_coeff(const double src[],
15                              SkScalar w,
16                              double coeff[3]) {
17    const double P20 = src[4] - src[0];
18    const double P10 = src[2] - src[0];
19    const double wP10 = w * P10;
20    coeff[0] = w * P20 - P20;
21    coeff[1] = P20 - 2 * wP10;
22    coeff[2] = wP10;
23}
24
25static double conic_eval_tan(const double coord[], SkScalar w, double t) {
26    double coeff[3];
27    conic_deriv_coeff(coord, w, coeff);
28    return t * (t * coeff[0] + coeff[1]) + coeff[2];
29}
30
31int SkDConic::FindExtrema(const double src[], SkScalar w, double t[1]) {
32    double coeff[3];
33    conic_deriv_coeff(src, w, coeff);
34
35    double tValues[2];
36    int roots = SkDQuad::RootsValidT(coeff[0], coeff[1], coeff[2], tValues);
37    // In extreme cases, the number of roots returned can be 2. Pathops
38    // will fail later on, so there's no advantage to plumbing in an error
39    // return here.
40    // SkASSERT(0 == roots || 1 == roots);
41
42    if (1 == roots) {
43        t[0] = tValues[0];
44        return 1;
45    }
46    return 0;
47}
48
49SkDVector SkDConic::dxdyAtT(double t) const {
50    SkDVector result = {
51        conic_eval_tan(&fPts[0].fX, fWeight, t),
52        conic_eval_tan(&fPts[0].fY, fWeight, t)
53    };
54    if (result.fX == 0 && result.fY == 0) {
55        if (zero_or_one(t)) {
56            result = fPts[2] - fPts[0];
57        } else {
58            // incomplete
59            SkDebugf("!k");
60        }
61    }
62    return result;
63}
64
65static double conic_eval_numerator(const double src[], SkScalar w, double t) {
66    SkASSERT(src);
67    SkASSERT(t >= 0 && t <= 1);
68    double src2w = src[2] * w;
69    double C = src[0];
70    double A = src[4] - 2 * src2w + C;
71    double B = 2 * (src2w - C);
72    return (A * t + B) * t + C;
73}
74
75
76static double conic_eval_denominator(SkScalar w, double t) {
77    double B = 2 * (w - 1);
78    double C = 1;
79    double A = -B;
80    return (A * t + B) * t + C;
81}
82
83bool SkDConic::hullIntersects(const SkDCubic& cubic, bool* isLinear) const {
84    return cubic.hullIntersects(*this, isLinear);
85}
86
87SkDPoint SkDConic::ptAtT(double t) const {
88    if (t == 0) {
89        return fPts[0];
90    }
91    if (t == 1) {
92        return fPts[2];
93    }
94    double denominator = conic_eval_denominator(fWeight, t);
95    SkDPoint result = {
96        conic_eval_numerator(&fPts[0].fX, fWeight, t) / denominator,
97        conic_eval_numerator(&fPts[0].fY, fWeight, t) / denominator
98    };
99    return result;
100}
101
102/* see quad subdivide for point rationale */
103/* w rationale : the mid point between t1 and t2 could be determined from the computed a/b/c
104   values if the computed w was known. Since we know the mid point at (t1+t2)/2, we'll assume
105   that it is the same as the point on the new curve t==(0+1)/2.
106
107    d / dz == conic_poly(dst, unknownW, .5) / conic_weight(unknownW, .5);
108
109    conic_poly(dst, unknownW, .5)
110                  =   a / 4 + (b * unknownW) / 2 + c / 4
111                  =  (a + c) / 4 + (bx * unknownW) / 2
112
113    conic_weight(unknownW, .5)
114                  =   unknownW / 2 + 1 / 2
115
116    d / dz                  == ((a + c) / 2 + b * unknownW) / (unknownW + 1)
117    d / dz * (unknownW + 1) ==  (a + c) / 2 + b * unknownW
118              unknownW       = ((a + c) / 2 - d / dz) / (d / dz - b)
119
120    Thus, w is the ratio of the distance from the mid of end points to the on-curve point, and the
121    distance of the on-curve point to the control point.
122 */
123SkDConic SkDConic::subDivide(double t1, double t2) const {
124    double ax, ay, az;
125    if (t1 == 0) {
126        ax = fPts[0].fX;
127        ay = fPts[0].fY;
128        az = 1;
129    } else if (t1 != 1) {
130        ax = conic_eval_numerator(&fPts[0].fX, fWeight, t1);
131        ay = conic_eval_numerator(&fPts[0].fY, fWeight, t1);
132        az = conic_eval_denominator(fWeight, t1);
133    } else {
134        ax = fPts[2].fX;
135        ay = fPts[2].fY;
136        az = 1;
137    }
138    double midT = (t1 + t2) / 2;
139    double dx = conic_eval_numerator(&fPts[0].fX, fWeight, midT);
140    double dy = conic_eval_numerator(&fPts[0].fY, fWeight, midT);
141    double dz = conic_eval_denominator(fWeight, midT);
142    double cx, cy, cz;
143    if (t2 == 1) {
144        cx = fPts[2].fX;
145        cy = fPts[2].fY;
146        cz = 1;
147    } else if (t2 != 0) {
148        cx = conic_eval_numerator(&fPts[0].fX, fWeight, t2);
149        cy = conic_eval_numerator(&fPts[0].fY, fWeight, t2);
150        cz = conic_eval_denominator(fWeight, t2);
151    } else {
152        cx = fPts[0].fX;
153        cy = fPts[0].fY;
154        cz = 1;
155    }
156    double bx = 2 * dx - (ax + cx) / 2;
157    double by = 2 * dy - (ay + cy) / 2;
158    double bz = 2 * dz - (az + cz) / 2;
159    SkDConic dst = {{{{ax / az, ay / az}, {bx / bz, by / bz}, {cx / cz, cy / cz}}
160            SkDEBUGPARAMS(fPts.fDebugGlobalState) },
161            SkDoubleToScalar(bz / sqrt(az * cz)) };
162    return dst;
163}
164
165SkDPoint SkDConic::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2,
166        SkScalar* weight) const {
167    SkDConic chopped = this->subDivide(t1, t2);
168    *weight = chopped.fWeight;
169    return chopped[1];
170}
171