1/*
2 * Copyright 2012 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 "CubicUtilities.h"
8#include "IntersectionUtilities.h"
9
10/*
11 Given a cubic c, t1, and t2, find a small cubic segment.
12
13 The new cubic is defined as points A, B, C, and D, where
14 s1 = 1 - t1
15 s2 = 1 - t2
16 A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1
17 D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2
18
19 We don't have B or C. So We define two equations to isolate them.
20 First, compute two reference T values 1/3 and 2/3 from t1 to t2:
21
22 c(at (2*t1 + t2)/3) == E
23 c(at (t1 + 2*t2)/3) == F
24
25 Next, compute where those values must be if we know the values of B and C:
26
27 _12   =  A*2/3 + B*1/3
28 12_   =  A*1/3 + B*2/3
29 _23   =  B*2/3 + C*1/3
30 23_   =  B*1/3 + C*2/3
31 _34   =  C*2/3 + D*1/3
32 34_   =  C*1/3 + D*2/3
33 _123  = (A*2/3 + B*1/3)*2/3 + (B*2/3 + C*1/3)*1/3 = A*4/9 + B*4/9 + C*1/9
34 123_  = (A*1/3 + B*2/3)*1/3 + (B*1/3 + C*2/3)*2/3 = A*1/9 + B*4/9 + C*4/9
35 _234  = (B*2/3 + C*1/3)*2/3 + (C*2/3 + D*1/3)*1/3 = B*4/9 + C*4/9 + D*1/9
36 234_  = (B*1/3 + C*2/3)*1/3 + (C*1/3 + D*2/3)*2/3 = B*1/9 + C*4/9 + D*4/9
37 _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3
38       =  A*8/27 + B*12/27 + C*6/27 + D*1/27
39       =  E
40 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3
41       =  A*1/27 + B*6/27 + C*12/27 + D*8/27
42       =  F
43 E*27  =  A*8    + B*12   + C*6     + D
44 F*27  =  A      + B*6    + C*12    + D*8
45
46Group the known values on one side:
47
48 M       = E*27 - A*8 - D     = B*12 + C* 6
49 N       = F*27 - A   - D*8   = B* 6 + C*12
50 M*2 - N = B*18
51 N*2 - M = C*18
52 B       = (M*2 - N)/18
53 C       = (N*2 - M)/18
54 */
55
56static double interp_cubic_coords(const double* src, double t)
57{
58    double ab = interp(src[0], src[2], t);
59    double bc = interp(src[2], src[4], t);
60    double cd = interp(src[4], src[6], t);
61    double abc = interp(ab, bc, t);
62    double bcd = interp(bc, cd, t);
63    double abcd = interp(abc, bcd, t);
64    return abcd;
65}
66
67void sub_divide(const Cubic& src, double t1, double t2, Cubic& dst) {
68    if (t1 == 0 && t2 == 1) {
69        dst[0] = src[0];
70        dst[1] = src[1];
71        dst[2] = src[2];
72        dst[3] = src[3];
73        return;
74    }
75    double ax = dst[0].x = interp_cubic_coords(&src[0].x, t1);
76    double ay = dst[0].y = interp_cubic_coords(&src[0].y, t1);
77    double ex = interp_cubic_coords(&src[0].x, (t1*2+t2)/3);
78    double ey = interp_cubic_coords(&src[0].y, (t1*2+t2)/3);
79    double fx = interp_cubic_coords(&src[0].x, (t1+t2*2)/3);
80    double fy = interp_cubic_coords(&src[0].y, (t1+t2*2)/3);
81    double dx = dst[3].x = interp_cubic_coords(&src[0].x, t2);
82    double dy = dst[3].y = interp_cubic_coords(&src[0].y, t2);
83    double mx = ex * 27 - ax * 8 - dx;
84    double my = ey * 27 - ay * 8 - dy;
85    double nx = fx * 27 - ax - dx * 8;
86    double ny = fy * 27 - ay - dy * 8;
87    /* bx = */ dst[1].x = (mx * 2 - nx) / 18;
88    /* by = */ dst[1].y = (my * 2 - ny) / 18;
89    /* cx = */ dst[2].x = (nx * 2 - mx) / 18;
90    /* cy = */ dst[2].y = (ny * 2 - my) / 18;
91}
92
93void sub_divide(const Cubic& src, const _Point& a, const _Point& d,
94        double t1, double t2, _Point dst[2]) {
95    double ex = interp_cubic_coords(&src[0].x, (t1 * 2 + t2) / 3);
96    double ey = interp_cubic_coords(&src[0].y, (t1 * 2 + t2) / 3);
97    double fx = interp_cubic_coords(&src[0].x, (t1 + t2 * 2) / 3);
98    double fy = interp_cubic_coords(&src[0].y, (t1 + t2 * 2) / 3);
99    double mx = ex * 27 - a.x * 8 - d.x;
100    double my = ey * 27 - a.y * 8 - d.y;
101    double nx = fx * 27 - a.x - d.x * 8;
102    double ny = fy * 27 - a.y - d.y * 8;
103    /* bx = */ dst[0].x = (mx * 2 - nx) / 18;
104    /* by = */ dst[0].y = (my * 2 - ny) / 18;
105    /* cx = */ dst[1].x = (nx * 2 - mx) / 18;
106    /* cy = */ dst[1].y = (ny * 2 - my) / 18;
107}
108
109/* classic one t subdivision */
110static void interp_cubic_coords(const double* src, double* dst, double t)
111{
112    double ab = interp(src[0], src[2], t);
113    double bc = interp(src[2], src[4], t);
114    double cd = interp(src[4], src[6], t);
115    double abc = interp(ab, bc, t);
116    double bcd = interp(bc, cd, t);
117    double abcd = interp(abc, bcd, t);
118
119    dst[0] = src[0];
120    dst[2] = ab;
121    dst[4] = abc;
122    dst[6] = abcd;
123    dst[8] = bcd;
124    dst[10] = cd;
125    dst[12] = src[6];
126}
127
128void chop_at(const Cubic& src, CubicPair& dst, double t)
129{
130    if (t == 0.5) {
131        dst.pts[0] = src[0];
132        dst.pts[1].x = (src[0].x + src[1].x) / 2;
133        dst.pts[1].y = (src[0].y + src[1].y) / 2;
134        dst.pts[2].x = (src[0].x + 2 * src[1].x + src[2].x) / 4;
135        dst.pts[2].y = (src[0].y + 2 * src[1].y + src[2].y) / 4;
136        dst.pts[3].x = (src[0].x + 3 * (src[1].x + src[2].x) + src[3].x) / 8;
137        dst.pts[3].y = (src[0].y + 3 * (src[1].y + src[2].y) + src[3].y) / 8;
138        dst.pts[4].x = (src[1].x + 2 * src[2].x + src[3].x) / 4;
139        dst.pts[4].y = (src[1].y + 2 * src[2].y + src[3].y) / 4;
140        dst.pts[5].x = (src[2].x + src[3].x) / 2;
141        dst.pts[5].y = (src[2].y + src[3].y) / 2;
142        dst.pts[6] = src[3];
143        return;
144    }
145    interp_cubic_coords(&src[0].x, &dst.pts[0].x, t);
146    interp_cubic_coords(&src[0].y, &dst.pts[0].y, t);
147}
148