OLD | NEW |
| (Empty) |
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 | |
46 Group 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 | |
56 static 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 | |
67 void 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 | |
93 void 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 */ | |
110 static 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 | |
128 void 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 } | |
OLD | NEW |