1 /*
2 * Copyright (c) 2007, 2017, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
23 * questions.
24 */
25
26 package sun.java2d.marlin;
27
28 final class DCurve {
29
30 double ax, ay, bx, by, cx, cy, dx, dy;
31 double dax, day, dbx, dby;
32
33 DCurve() {
34 }
35
36 void set(double[] points, int type) {
37 switch(type) {
38 case 8:
39 set(points[0], points[1],
40 points[2], points[3],
41 points[4], points[5],
42 points[6], points[7]);
43 return;
44 case 6:
45 set(points[0], points[1],
46 points[2], points[3],
47 points[4], points[5]);
48 return;
49 default:
50 throw new InternalError("Curves can only be cubic or quadratic");
51 }
52 }
53
54 void set(double x1, double y1,
55 double x2, double y2,
56 double x3, double y3,
57 double x4, double y4)
58 {
59 ax = 3.0d * (x2 - x3) + x4 - x1;
60 ay = 3.0d * (y2 - y3) + y4 - y1;
61 bx = 3.0d * (x1 - 2.0d * x2 + x3);
62 by = 3.0d * (y1 - 2.0d * y2 + y3);
63 cx = 3.0d * (x2 - x1);
64 cy = 3.0d * (y2 - y1);
65 dx = x1;
66 dy = y1;
67 dax = 3.0d * ax; day = 3.0d * ay;
68 dbx = 2.0d * bx; dby = 2.0d * by;
69 }
70
71 void set(double x1, double y1,
72 double x2, double y2,
73 double x3, double y3)
74 {
75 ax = 0.0d; ay = 0.0d;
76 bx = x1 - 2.0d * x2 + x3;
77 by = y1 - 2.0d * y2 + y3;
78 cx = 2.0d * (x2 - x1);
79 cy = 2.0d * (y2 - y1);
80 dx = x1;
81 dy = y1;
82 dax = 0.0d; day = 0.0d;
83 dbx = 2.0d * bx; dby = 2.0d * by;
84 }
85
86 double xat(double t) {
87 return t * (t * (t * ax + bx) + cx) + dx;
88 }
89 double yat(double t) {
90 return t * (t * (t * ay + by) + cy) + dy;
91 }
92
93 double dxat(double t) {
94 return t * (t * dax + dbx) + cx;
95 }
96
97 double dyat(double t) {
98 return t * (t * day + dby) + cy;
99 }
100
101 int dxRoots(double[] roots, int off) {
102 return DHelpers.quadraticRoots(dax, dbx, cx, roots, off);
103 }
104
105 int dyRoots(double[] roots, int off) {
106 return DHelpers.quadraticRoots(day, dby, cy, roots, off);
107 }
108
109 int infPoints(double[] pts, int off) {
110 // inflection point at t if -f'(t)x*f''(t)y + f'(t)y*f''(t)x == 0
111 // Fortunately, this turns out to be quadratic, so there are at
112 // most 2 inflection points.
113 final double a = dax * dby - dbx * day;
114 final double b = 2.0d * (cy * dax - day * cx);
115 final double c = cy * dbx - cx * dby;
116
117 return DHelpers.quadraticRoots(a, b, c, pts, off);
118 }
119
120 // finds points where the first and second derivative are
121 // perpendicular. This happens when g(t) = f'(t)*f''(t) == 0 (where
122 // * is a dot product). Unfortunately, we have to solve a cubic.
123 private int perpendiculardfddf(double[] pts, int off) {
124 assert pts.length >= off + 4;
125
126 // these are the coefficients of some multiple of g(t) (not g(t),
127 // because the roots of a polynomial are not changed after multiplication
128 // by a constant, and this way we save a few multiplications).
129 final double a = 2.0d * (dax*dax + day*day);
130 final double b = 3.0d * (dax*dbx + day*dby);
131 final double c = 2.0d * (dax*cx + day*cy) + dbx*dbx + dby*dby;
132 final double d = dbx*cx + dby*cy;
133 return DHelpers.cubicRootsInAB(a, b, c, d, pts, off, 0.0d, 1.0d);
134 }
135
136 // Tries to find the roots of the function ROC(t)-w in [0, 1). It uses
137 // a variant of the false position algorithm to find the roots. False
138 // position requires that 2 initial values x0,x1 be given, and that the
139 // function must have opposite signs at those values. To find such
140 // values, we need the local extrema of the ROC function, for which we
141 // need the roots of its derivative; however, it's harder to find the
142 // roots of the derivative in this case than it is to find the roots
143 // of the original function. So, we find all points where this curve's
144 // first and second derivative are perpendicular, and we pretend these
145 // are our local extrema. There are at most 3 of these, so we will check
146 // at most 4 sub-intervals of (0,1). ROC has asymptotes at inflection
147 // points, so roc-w can have at least 6 roots. This shouldn't be a
148 // problem for what we're trying to do (draw a nice looking curve).
149 int rootsOfROCMinusW(double[] roots, int off, final double w, final double err) {
150 // no OOB exception, because by now off<=6, and roots.length >= 10
151 assert off <= 6 && roots.length >= 10;
152 int ret = off;
153 int numPerpdfddf = perpendiculardfddf(roots, off);
154 double t0 = 0.0d, ft0 = ROCsq(t0) - w*w;
155 roots[off + numPerpdfddf] = 1.0d; // always check interval end points
156 numPerpdfddf++;
157 for (int i = off; i < off + numPerpdfddf; i++) {
158 double t1 = roots[i], ft1 = ROCsq(t1) - w*w;
159 if (ft0 == 0.0d) {
160 roots[ret++] = t0;
161 } else if (ft1 * ft0 < 0.0d) { // have opposite signs
162 // (ROC(t)^2 == w^2) == (ROC(t) == w) is true because
163 // ROC(t) >= 0 for all t.
164 roots[ret++] = falsePositionROCsqMinusX(t0, t1, w*w, err);
165 }
166 t0 = t1;
167 ft0 = ft1;
168 }
169
170 return ret - off;
171 }
172
173 private static double eliminateInf(double x) {
174 return (x == Double.POSITIVE_INFINITY ? Double.MAX_VALUE :
175 (x == Double.NEGATIVE_INFINITY ? Double.MIN_VALUE : x));
176 }
177
178 // A slight modification of the false position algorithm on wikipedia.
179 // This only works for the ROCsq-x functions. It might be nice to have
180 // the function as an argument, but that would be awkward in java6.
181 // TODO: It is something to consider for java8 (or whenever lambda
182 // expressions make it into the language), depending on how closures
183 // and turn out. Same goes for the newton's method
184 // algorithm in DHelpers.java
185 private double falsePositionROCsqMinusX(double x0, double x1,
186 final double x, final double err)
187 {
188 final int iterLimit = 100;
189 int side = 0;
190 double t = x1, ft = eliminateInf(ROCsq(t) - x);
191 double s = x0, fs = eliminateInf(ROCsq(s) - x);
192 double r = s, fr;
193 for (int i = 0; i < iterLimit && Math.abs(t - s) > err * Math.abs(t + s); i++) {
194 r = (fs * t - ft * s) / (fs - ft);
195 fr = ROCsq(r) - x;
196 if (sameSign(fr, ft)) {
197 ft = fr; t = r;
198 if (side < 0) {
199 fs /= (1 << (-side));
200 side--;
201 } else {
202 side = -1;
203 }
204 } else if (fr * fs > 0) {
205 fs = fr; s = r;
206 if (side > 0) {
207 ft /= (1 << side);
208 side++;
209 } else {
210 side = 1;
211 }
212 } else {
213 break;
214 }
215 }
216 return r;
217 }
218
219 private static boolean sameSign(double x, double y) {
220 // another way is to test if x*y > 0. This is bad for small x, y.
221 return (x < 0.0d && y < 0.0d) || (x > 0.0d && y > 0.0d);
222 }
223
224 // returns the radius of curvature squared at t of this curve
225 // see http://en.wikipedia.org/wiki/Radius_of_curvature_(applications)
226 private double ROCsq(final double t) {
227 // dx=xat(t) and dy=yat(t). These calls have been inlined for efficiency
228 final double dx = t * (t * dax + dbx) + cx;
229 final double dy = t * (t * day + dby) + cy;
230 final double ddx = 2.0d * dax * t + dbx;
231 final double ddy = 2.0d * day * t + dby;
232 final double dx2dy2 = dx*dx + dy*dy;
233 final double ddx2ddy2 = ddx*ddx + ddy*ddy;
234 final double ddxdxddydy = ddx*dx + ddy*dy;
235 return dx2dy2*((dx2dy2*dx2dy2) / (dx2dy2 * ddx2ddy2 - ddxdxddydy*ddxdxddydy));
236 }
237 }
--- EOF ---