1 /*
2 * jcurl java curling software framework http://www.jcurl.orgCopyright (C)
3 * 2005-2009 M. Rohrmoser
4 *
5 * This program is free software; you can redistribute it and/or modify it under
6 * the terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 2 of the License, or (at your option) any later
8 * version.
9 *
10 * This program is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License along with
16 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple
17 * Place, Suite 330, Boston, MA 02111-1307 USA
18 */
19 package org.jcurl.math;
20
21 import java.awt.Shape;
22 import java.awt.geom.GeneralPath;
23
24 import org.apache.commons.logging.Log;
25 import org.jcurl.core.impl.CurlerDenny;
26 import org.jcurl.core.log.JCLoggerFactory;
27
28 /**
29 * Helper for convenient approximated Java2D drawing of arbitratry
30 * {@link R1RNFunction}s with at least 2 dimensions.
31 *
32 * @see Shaper
33 * @see Shapeable
34 * @author <a href="mailto:m@jcurl.org">M. Rohrmoser </a>
35 * @version $Id: ShaperUtils.java 1031 2009-07-23 15:06:05Z mro $
36 */
37 public abstract class ShaperUtils {
38
39 private static final Log log = JCLoggerFactory.getLogger(ShaperUtils.class);
40
41 /**
42 * Compute the control points and add one <a
43 * href="http://en.wikipedia.org/wiki/B%C3%A9zier_curve">Cubic Bezier Curve</a>
44 * to a {@link GeneralPath}. Does <b>no</b> initial
45 * {@link GeneralPath#moveTo(float, float)}.
46 *
47 * <h3>Approximation algorithm</h3>
48 * <p>
49 * This ansatz uses no adaptive optimisation but the nature of curves as
50 * they're typical to curling:
51 * <ul>
52 * <li>continuous within [tmin:tmax] - at least C0, C1</li>
53 * <li>smoothly increasing curvature</li>
54 * <li>not meandering</li>
55 * </ul>
56 * So we use
57 * <ul>
58 * <li>the start- and endpoint of each interval als control points k0 and
59 * k3</li>
60 * <li>the directions (normalised velocities) in the control points k0 and
61 * k3</li>
62 * <li>a ratio 3:2:1 of the distances |k0-k1| : |k1-k2| : |k2-k3|</li>
63 * </ul>
64 * <p>
65 * This causes quite a computation - without iteration/recursion though, but
66 * 1 square root and many double multiplications - but this is well worth
67 * while as we can reduce the curve segments to draw significantly. One
68 * cubic bezier curve per seven meters curve length gives an error < 2 mm
69 * (using {@link CurlerDenny} with 24s draw-to-tee and 1m curl)!
70 * </p>
71 * <p>
72 * TODO maybe re-use endpoint location and velocity. This can cause pain at
73 * C1 discontinuous t's (collissions).
74 * </p>
75 * <h3><a href="http://en.wikipedia.org/wiki/Maxima_(software)">Maxima</a>
76 * Solution</h3>
77 *
78 * <pre>
79 * radsubstflag: true$
80 * k1_0 = k0_0 + l * v0_0;
81 * k1_1 = k0_1 + l * v0_1;
82 * k2_0 = k3_0 - n * v3_0;
83 * k2_1 = k3_1 - n * v3_1;
84 * l/n=a/c;
85 * ((k2_0 - k1_0)*(k2_0 - k1_0) + (k2_1 - k1_1)*(k2_1 - k1_1)) / (n*n) = b*b / (c*c);
86 * solve([%th(6), %th(5), %th(4), %th(3), %th(2), %th(1)],[k1_0, k1_1, k2_0, k2_1, l, n]);
87 * factor(%);
88 * ratsimp(%);
89 * ratsubst(V0, v0_1ˆ2+v0_0ˆ2, %);
90 * ratsubst(V3, v3_1ˆ2+v3_0ˆ2, %);
91 * ratsubst(A, k0_1-k3_1, %);
92 * ratsubst(B, k0_0-k3_0, %);
93 * ratsubst(T, 2*a*c*v0_0*v3_0+aˆ2*v0_1ˆ2+aˆ2*v0_0ˆ2-bˆ2, %);
94 * ratsubst(Q, cˆ2*V3+aˆ2*V0+T+2*a*c*v0_1*v3_1-aˆ2*v0_1ˆ2-aˆ2*v0_0ˆ2, %);
95 * ratsubst(W, Bˆ2*T+Bˆ2*(bˆ2-Q)+cˆ2*(v3_0ˆ2*Bˆ2-v3_0ˆ2*Aˆ2)-aˆ2*v0_1ˆ2*Bˆ2+v3_1*(2*cˆ2*v3_0*A*B
96 * +2*a*c*v0_0*A*B)+v0_1*(2*a*c*v3_0*A*B+2*aˆ2*v0_0*A*B)-2*a*c*v0_0*v3_0*Aˆ2-aˆ2*v0_0ˆ2*Aˆ2
97 * +bˆ2*Aˆ2, %);
98 * expand(%);
99 * factor(%);
100 * ratsubst(R, c*v3_0*B+a*v0_0*B+c*v3_1*A+a*v0_1*A, %);
101 * </pre>
102 */
103 static void curveTo(final R1RNFunction f, final double tmin,
104 final double tmax, final GeneralPath gp, final float zoom) {
105 final double eps = 1e-6;
106
107 // first control point (startpoint). The same as gp.getCurrentPoint()
108 final double k0_0 = f.at(tmin, 0, 0);
109 final double k0_1 = f.at(tmin, 0, 1);
110 // normalized startpoint velocity
111 double v0_0 = f.at(tmin, 1, 0);
112 double v0_1 = f.at(tmin, 1, 1);
113 if (v0_0 * v0_0 + v0_1 * v0_1 < eps) {
114 v0_0 = f.at(tmin + eps, 1, 0);
115 v0_1 = f.at(tmin + eps, 1, 1);
116 }
117 double v = Math.sqrt(v0_0 * v0_0 + v0_1 * v0_1);
118 v0_0 /= v;
119 v0_1 /= v;
120
121 // 4th control point (endpoint).
122 final double k3_0 = f.at(tmax, 0, 0);
123 final double k3_1 = f.at(tmax, 0, 1);
124 // normalized endpoint velocity
125 double v3_0 = f.at(tmax, 1, 0);
126 double v3_1 = f.at(tmax, 1, 1);
127 if (v3_0 * v3_0 + v3_1 * v3_1 < eps) {
128 v3_0 = f.at(tmax - eps, 1, 0);
129 v3_1 = f.at(tmax - eps, 1, 1);
130 }
131 v = Math.sqrt(v3_0 * v3_0 + v3_1 * v3_1);
132 v3_0 /= v;
133 v3_1 /= v;
134
135 final double a = 3;
136 final double b = 2;
137 final double c = 1;
138 final double V0 = v0_1 * v0_1 + v0_0 * v0_0;
139 final double V3 = v3_1 * v3_1 + v3_0 * v3_0;
140 final double A = k0_1 - k3_1;
141 final double B = k0_0 - k3_0;
142 final double T = 2 * a * c * v0_0 * v3_0 + a * a * v0_1 * v0_1 + a * a
143 * v0_0 * v0_0 - b * b;
144 final double Q = c * c * V3 + a * a * V0 + T + 2 * a * c * v0_1 * v3_1
145 - a * a * v0_1 * v0_1 - a * a * v0_0 * v0_0;
146 double W = B * B * T + B * B * (b * b - Q) + c * c
147 * (v3_0 * v3_0 * B * B - v3_0 * v3_0 * A * A) - a * a * v0_1
148 * v0_1 * B * B + v3_1 * 2 * c * c * v3_0 * A * B + 2 * a * c
149 * v0_0 * A * B + v0_1
150 * (2 * a * c * v3_0 * A * B + 2 * a * a * v0_0 * A * B) - 2 * a
151 * c * v0_0 * v3_0 * A * A - a * a * v0_0 * v0_0 * A * A + b * b
152 * A * A;
153 if (W < 0) {
154 if (log.isWarnEnabled()) {
155 log.warn("Arithmetic trouble:");
156 log.warn("v0=(" + v0_0 + ", " + v0_1 + ")");
157 log.warn("v3=(" + v3_0 + ", " + v3_1 + ")");
158 log.warn("V0=" + V0);
159 log.warn("V3=" + V3);
160 log.warn("A=" + A);
161 log.warn("B=" + B);
162 log.warn("T=" + T);
163 log.warn("Q=" + Q);
164 log.warn("W=" + W);
165 }
166 gp.moveTo(zoom * (float) k3_0, zoom * (float) k3_1);
167 return;
168 }
169 W = Math.sqrt(W);
170 final double R = c * v3_0 * B + a * v0_0 * B + c * v3_1 * A + a * v0_1
171 * A;
172
173 final double l, n;
174 if (true) {
175 final double F = (W + R) / Q;
176 l = -a * F;
177 n = -c * F;
178 } else {
179 final double F = (W - R) / Q;
180 l = a * F;
181 n = c * F;
182 }
183 if (Double.isNaN(l) || Double.isNaN(n)) {
184 log.warn("v0=(" + v0_0 + ", " + v0_1 + ")");
185 log.warn("v3=(" + v3_0 + ", " + v3_1 + ")");
186 log.warn("V0=" + V0);
187 log.warn("V3=" + V3);
188 log.warn("A=" + A);
189 log.warn("B=" + B);
190 log.warn("T=" + T);
191 log.warn("Q=" + Q);
192 log.warn("W=" + W);
193 log.warn("R=" + R);
194 }
195
196 final float k1_0 = (float) (k0_0 + l * v0_0);
197 final float k1_1 = (float) (k0_1 + l * v0_1);
198 final float k2_0 = (float) (k3_0 - n * v3_0);
199 final float k2_1 = (float) (k3_1 - n * v3_1);
200 if (log.isDebugEnabled())
201 log.debug("(" + k1_0 + ", " + k1_1 + "), (" + k2_0 + ", " + k2_1
202 + "), (" + (float) k3_0 + ", " + (float) k3_1 + ")");
203 gp.curveTo(zoom * k1_0, zoom * k1_1, zoom * k2_0, zoom * k2_1, zoom
204 * (float) k3_0, zoom * (float) k3_1);
205 }
206
207 private static float interpolate(final float min, final float max,
208 final float t, final Interpolator ip) {
209 final float d = max - min;
210 return min + d * ip.interpolate((t - min) / d);
211 }
212
213 /**
214 * Interpolate using <a
215 * href="http://en.wikipedia.org/wiki/B%C3%A9zier_curve">Cubic Bezier Curves</a>.
216 * <p>
217 * Computes the required intermediate <code>t</code> samples and delegates
218 * to {@link #curveTo(R1RNFunction, double, double, GeneralPath, float)} to
219 * compute the interpolating curve segments.
220 * </p>
221 *
222 * @param src
223 * the (at least 2-dimensional) curve. Higher dimensions are
224 * ignored.
225 * @param min
226 * the min input <code>t</code> to
227 * {@link R1RNFunction#at(double, int, int)}
228 * @param max
229 * the max input <code>t</code> to
230 * {@link R1RNFunction#at(double, int, int)}
231 * @param curves
232 * the number of interpolating cubic bezier curves - must be
233 * >= 1.
234 * @param zoom
235 * graphics zoom factor (typically 1)
236 * @param ip
237 * the {@link Interpolator} to get the intermediate
238 * <code>t</code> sample values.
239 * @see #curveTo(R1RNFunction, double, double, GeneralPath, float)
240 */
241 public static Shape interpolateCubic(final R1RNFunction src,
242 final double min, final double max, final int curves,
243 final float zoom, final Interpolator ip) {
244 // setup
245 if (curves < 1)
246 throw new IllegalArgumentException(
247 "Give me at least 1 (connect start + stop)");
248 final float d = (float) (max - min);
249 final GeneralPath gp = new GeneralPath(GeneralPath.WIND_NON_ZERO,
250 3 * curves + 1); // +1 just to be sure...
251 // start
252 final float x = (float) src.at(min, 0, 0);
253 final float y = (float) src.at(min, 0, 1);
254 gp.moveTo(zoom * x, zoom * y);
255
256 double told = min;
257 // intermediate
258 final int n = curves;
259 for (int i = 1; i < n; i++) {
260 final double t = min + d * ip.interpolate((float) i / n);
261 curveTo(src, told, t, gp, zoom);
262 told = t;
263 }
264
265 // stop
266 curveTo(src, told, max, gp, zoom);
267 return gp;
268 }
269
270 /**
271 * Interpolate using <a
272 * href="http://en.wikipedia.org/wiki/B%C3%A9zier_curve">Linear Bezier
273 * Curves</a>.
274 * <p>
275 * Computes the required intermediate <code>t</code> samples and delegates
276 * to {@link #lineTo(R1RNFunction, double, GeneralPath, float)} to compute
277 * the interpolating curve segments.
278 * </p>
279 *
280 * @param src
281 * the (at least 2-dimensional) curve. Higher dimensions are
282 * ignored.
283 * @param min
284 * the min input <code>t</code> to
285 * {@link R1RNFunction#at(double, int, int)}
286 * @param max
287 * the max input <code>t</code> to
288 * {@link R1RNFunction#at(double, int, int)}
289 * @param curves
290 * the number of line segments - must be >= 1.
291 * @param zoom
292 * graphics zoom factor (typically 1)
293 * @param ip
294 * the {@link Interpolator} to get the intermediate sample
295 * <code>t</code> values.
296 * @see #lineTo(R1RNFunction, double, GeneralPath, float)
297 */
298 public static Shape interpolateLinear(final R1RNFunction src,
299 final double min, final double max, final int curves,
300 final float zoom, final Interpolator ip) {
301 // setup
302 if (curves < 1)
303 throw new IllegalArgumentException(
304 "Give me at least 1 (connect start + stop)");
305 final float d = (float) (max - min);
306 final GeneralPath gp = new GeneralPath(GeneralPath.WIND_NON_ZERO,
307 curves + 1); // +1 just to be sure...
308 // start
309 final float x = (float) src.at(min, 0, 0);
310 final float y = (float) src.at(min, 0, 1);
311 gp.moveTo(zoom * x, zoom * y);
312
313 // intermediate
314 final int n = curves;
315 for (int i = 1; i < n; i++) {
316 final double t = min + d * ip.interpolate((float) i / n);
317 lineTo(src, t, gp, zoom);
318 }
319
320 // stop
321 lineTo(src, max, gp, zoom);
322 return gp;
323 }
324
325 /**
326 * Interpolate using <a
327 * href="http://en.wikipedia.org/wiki/B%C3%A9zier_curve">Quadratic Bezier
328 * Curves</a>.
329 * <p>
330 * Computes the required intermediate <code>t</code> samples and delegates
331 * to {@link #quadTo(R1RNFunction, double, double, GeneralPath, float)} to
332 * compute the interpolating curve segments.
333 * </p>
334 *
335 * @param src
336 * the (2-dimensional) curve. Higher dimensions are ignored.
337 * @param min
338 * the min input <code>t</code> to
339 * {@link R1RNFunction#at(double, int, int)}
340 * @param max
341 * the max input <code>t</code> to
342 * {@link R1RNFunction#at(double, int, int)}
343 * @param curves
344 * the number of line segments - must be >= 1.
345 * @param zoom
346 * graphics zoom factor (typically 1)
347 * @param ip
348 * the {@link Interpolator} to get the intermediate sample
349 * <code>t</code> values.
350 * @see #quadTo(R1RNFunction, double, double, GeneralPath, float)
351 */
352 public static Shape interpolateQuadratic(final R1RNFunction src,
353 final double min, final double max, final int curves,
354 final float zoom, final Interpolator ip) {
355 // setup
356 if (curves < 1)
357 throw new IllegalArgumentException(
358 "Give me at least 1 (connect start + stop)");
359 final float d = (float) (max - min);
360 final GeneralPath gp = new GeneralPath(GeneralPath.WIND_NON_ZERO,
361 2 * curves + 1); // +1 just to be sure...
362 // start
363 final float x = (float) src.at(min, 0, 0);
364 final float y = (float) src.at(min, 0, 1);
365 gp.moveTo(zoom * x, zoom * y);
366
367 double told = min;
368 // intermediate
369 final int n = curves;
370 for (int i = 1; i < n; i++) {
371 final double t = min + d * ip.interpolate((float) i / n);
372 quadTo(src, told, t, gp, zoom);
373 told = t;
374 }
375
376 // stop
377 quadTo(src, told, max, gp, zoom);
378 return gp;
379 }
380
381 /**
382 * Add one <a href="http://en.wikipedia.org/wiki/B%C3%A9zier_curve">Linear
383 * Bezier Curve</a> to a {@link GeneralPath}. Does <b>no</b> initial
384 * {@link GeneralPath#moveTo(float, float)}.
385 *
386 * <h3>Approximation algorithm</h3>
387 * <p>
388 * Just connect start- and endpoint.
389 * <p>
390 * TODO maybe re-use endpoint location and velocity. This can cause pain at
391 * C1 discontinuous t's (collissions).
392 * </p>
393 */
394 static final void lineTo(final R1RNFunction f, final double tmax,
395 final GeneralPath gp, final float zoom) {
396 final float x = (float) f.at(tmax, 0, 0);
397 final float y = (float) f.at(tmax, 0, 1);
398 gp.lineTo(zoom * x, zoom * y);
399 }
400
401 /**
402 * Compute the control point and add one <a
403 * href="http://en.wikipedia.org/wiki/B%C3%A9zier_curve">Quadratic Bezier
404 * Curve</a> to a {@link GeneralPath}. Does <b>no</b> initial
405 * {@link GeneralPath#moveTo(float, float)}.
406 *
407 * <h3>Approximation algorithm</h3>
408 * <p>
409 * This ansatz uses no adaptive optimisation but only
410 * <ul>
411 * <li>the start- and endpoint of each interval als control points k0 and
412 * k2</li>
413 * <li>the directions (normalised velocities) in the control points k0 and
414 * k2. The intersection is used as k1.</li>
415 * </ul>
416 * <p>
417 * TODO maybe re-use endpoint location and velocity. This can cause pain at
418 * C1 discontinuous t's (collissions).
419 * </p>
420 * <h3><a href="http://en.wikipedia.org/wiki/Maxima_(software)">Maxima</a>
421 * Solution</h3>
422 *
423 * <pre>
424 * radsubstflag: true$
425 * k0_0 + l * v0_0 = k2_0 + m * v2_0;
426 * k0_1 + l * v0_1 = k2_1 + m * v2_1;
427 * solve([%th(2),%th(1)],[l,m]);
428 * subst(q, v0_1 * v2_0 - v0_0 * v2_1, %);
429 * subst(dx_0 + k0_0, k2_0, %);
430 * subst(dx_1 + k0_1, k2_1, %);
431 * ratsimp(%);
432 * </pre>
433 */
434 static final void quadTo(final R1RNFunction f, final double tmin,
435 final double tmax, final GeneralPath gp, final float zoom) {
436 final double eps = 1e-6;
437
438 // first control point (startpoint). The same as gp.getCurrentPoint()
439 final double k0_0 = f.at(tmin, 0, 0);
440 final double k0_1 = f.at(tmin, 0, 1);
441 // startpoint velocity
442 double v0_0 = f.at(tmin, 1, 0);
443 double v0_1 = f.at(tmin, 1, 1);
444 if (v0_0 * v0_0 + v0_1 * v0_1 < eps) {
445 v0_0 = f.at(tmin + eps, 1, 0);
446 v0_1 = f.at(tmin + eps, 1, 1);
447 }
448
449 // 3rd control point (endpoint).
450 final double k2_0 = f.at(tmax, 0, 0);
451 final double k2_1 = f.at(tmax, 0, 1);
452 // endpoint velocity
453 double v2_0 = f.at(tmax, 1, 0);
454 double v2_1 = f.at(tmax, 1, 1);
455 if (v2_0 * v2_0 + v2_1 * v2_1 < eps) {
456 v2_0 = f.at(tmax - eps, 1, 0);
457 v2_1 = f.at(tmax - eps, 1, 1);
458 }
459
460 // compute the 2nd control point
461 final double dx_0 = k2_0 - k0_0;
462 final double dx_1 = k2_1 - k0_1;
463 final double q = v0_1 * v2_0 - v0_0 * v2_1;
464 final double m = -(dx_0 * v0_1 - dx_1 * v0_0) / q;
465
466 // 2nd control point is
467 final float k1_0 = (float) (k2_0 + m * v2_0);
468 final float k1_1 = (float) (k2_1 + m * v2_1);
469
470 if (true)
471 gp.quadTo(zoom * k1_0, zoom * k1_1, zoom * (float) k2_0, zoom
472 * (float) k2_1);
473 else {
474 gp.lineTo(zoom * k1_0, zoom * k1_1);
475 gp.lineTo(zoom * (float) k2_0, zoom * (float) k2_1);
476 }
477 }
478
479 static String toString(final double[] arr) {
480 final StringBuilder w = new StringBuilder();
481 if (arr == null)
482 w.append("null");
483 else {
484 boolean start = true;
485 w.append("[");
486 for (final double element : arr) {
487 if (!start)
488 w.append(" ");
489 w.append(Double.toString(element));
490 start = false;
491 }
492 w.append("]");
493 }
494 return w.toString();
495 }
496 }