001/* 002 * @(#)CubicCurve2D.java 1.35 06/04/17 003 * 004 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 005 * SUN PROPRIETARY/CONFIDENTIAL. Use is subject to license terms. 006 */ 007 008package armyc2.c2sd.graphics2d; 009 010/** 011 * The <code>CubicCurve2D</code> class defines a cubic parametric curve 012 * segment in {@code (x,y)} coordinate space. 013 * <p> 014 * This class is only the abstract superclass for all objects which 015 * store a 2D cubic curve segment. 016 * The actual storage representation of the coordinates is left to 017 * the subclass. 018 * 019 * @version 1.35, 04/17/06 020 * @author Jim Graham 021 * @since 1.2 022 */ 023public /* abstract */ final class CubicCurve2D /*implements Shape, Cloneable*/ { 024 025 /** 026 * This is an abstract class that cannot be instantiated directly. 027 * Type-specific implementation subclasses are available for 028 * instantiation and provide a number of formats for storing 029 * the information necessary to satisfy the various accessor 030 * methods below. 031 * 032 * @see java.awt.geom.CubicCurve2D.Float 033 * @see java.awt.geom.CubicCurve2D.Double 034 * @since 1.2 035 */ 036// protected CubicCurve2D() { 037// } 038 039 /** 040 * Returns the square of the flatness of the cubic curve specified 041 * by the indicated control points. The flatness is the maximum distance 042 * of a control point from the line connecting the end points. 043 * 044 * @param x1 the X coordinate that specifies the start point 045 * of a {@code CubicCurve2D} 046 * @param y1 the Y coordinate that specifies the start point 047 * of a {@code CubicCurve2D} 048 * @param ctrlx1 the X coordinate that specifies the first control point 049 * of a {@code CubicCurve2D} 050 * @param ctrly1 the Y coordinate that specifies the first control point 051 * of a {@code CubicCurve2D} 052 * @param ctrlx2 the X coordinate that specifies the second control point 053 * of a {@code CubicCurve2D} 054 * @param ctrly2 the Y coordinate that specifies the second control point 055 * of a {@code CubicCurve2D} 056 * @param x2 the X coordinate that specifies the end point 057 * of a {@code CubicCurve2D} 058 * @param y2 the Y coordinate that specifies the end point 059 * of a {@code CubicCurve2D} 060 * @return the square of the flatness of the {@code CubicCurve2D} 061 * represented by the specified coordinates. 062 * @since 1.2 063 */ 064 public static double getFlatnessSq2(double x1, double y1, 065 double ctrlx1, double ctrly1, 066 double ctrlx2, double ctrly2, 067 double x2, double y2) { 068 //return Math.max(Line2D.ptSegDistSq(x1, y1, x2, y2, ctrlx1, ctrly1), 069 // Line2D.ptSegDistSq(x1, y1, x2, y2, ctrlx2, ctrly2)); 070 071 return Math.max(Line2D.ptLineDistSq(x1, y1, x2, y2, ctrlx1, ctrly1), 072 Line2D.ptLineDistSq(x1, y1, x2, y2, ctrlx2, ctrly2)); 073 } 074 075 /** 076 * Returns the flatness of the cubic curve specified 077 * by the indicated control points. The flatness is the maximum distance 078 * of a control point from the line connecting the end points. 079 * 080 * @param x1 the X coordinate that specifies the start point 081 * of a {@code CubicCurve2D} 082 * @param y1 the Y coordinate that specifies the start point 083 * of a {@code CubicCurve2D} 084 * @param ctrlx1 the X coordinate that specifies the first control point 085 * of a {@code CubicCurve2D} 086 * @param ctrly1 the Y coordinate that specifies the first control point 087 * of a {@code CubicCurve2D} 088 * @param ctrlx2 the X coordinate that specifies the second control point 089 * of a {@code CubicCurve2D} 090 * @param ctrly2 the Y coordinate that specifies the second control point 091 * of a {@code CubicCurve2D} 092 * @param x2 the X coordinate that specifies the end point 093 * of a {@code CubicCurve2D} 094 * @param y2 the Y coordinate that specifies the end point 095 * of a {@code CubicCurve2D} 096 * @return the flatness of the {@code CubicCurve2D} 097 * represented by the specified coordinates. 098 * @since 1.2 099 */ 100 public static double getFlatness(double x1, double y1, 101 double ctrlx1, double ctrly1, 102 double ctrlx2, double ctrly2, 103 double x2, double y2) { 104 return Math.sqrt(getFlatnessSq2(x1, y1, ctrlx1, ctrly1, 105 ctrlx2, ctrly2, x2, y2)); 106 } 107 108 /** 109 * Returns the square of the flatness of the cubic curve specified 110 * by the control points stored in the indicated array at the 111 * indicated index. The flatness is the maximum distance 112 * of a control point from the line connecting the end points. 113 * @param coords an array containing coordinates 114 * @param offset the index of <code>coords</code> from which to begin 115 * getting the end points and control points of the curve 116 * @return the square of the flatness of the <code>CubicCurve2D</code> 117 * specified by the coordinates in <code>coords</code> at 118 * the specified offset. 119 * @since 1.2 120 */ 121 public static double getFlatnessSq(double coords[], int offset) { 122 return getFlatnessSq2(coords[offset + 0], coords[offset + 1], 123 coords[offset + 2], coords[offset + 3], 124 coords[offset + 4], coords[offset + 5], 125 coords[offset + 6], coords[offset + 7]); 126 } 127 128 /** 129 * Returns the flatness of the cubic curve specified 130 * by the control points stored in the indicated array at the 131 * indicated index. The flatness is the maximum distance 132 * of a control point from the line connecting the end points. 133 * @param coords an array containing coordinates 134 * @param offset the index of <code>coords</code> from which to begin 135 * getting the end points and control points of the curve 136 * @return the flatness of the <code>CubicCurve2D</code> 137 * specified by the coordinates in <code>coords</code> at 138 * the specified offset. 139 * @since 1.2 140 */ 141 public static double getFlatness2(double coords[], int offset) { 142 return getFlatness(coords[offset + 0], coords[offset + 1], 143 coords[offset + 2], coords[offset + 3], 144 coords[offset + 4], coords[offset + 5], 145 coords[offset + 6], coords[offset + 7]); 146 } 147 148 /** 149 * Subdivides the cubic curve specified by the coordinates 150 * stored in the <code>src</code> array at indices <code>srcoff</code> 151 * through (<code>srcoff</code> + 7) and stores the 152 * resulting two subdivided curves into the two result arrays at the 153 * corresponding indices. 154 * Either or both of the <code>left</code> and <code>right</code> 155 * arrays may be <code>null</code> or a reference to the same array 156 * as the <code>src</code> array. 157 * Note that the last point in the first subdivided curve is the 158 * same as the first point in the second subdivided curve. Thus, 159 * it is possible to pass the same array for <code>left</code> 160 * and <code>right</code> and to use offsets, such as <code>rightoff</code> 161 * equals (<code>leftoff</code> + 6), in order 162 * to avoid allocating extra storage for this common point. 163 * @param src the array holding the coordinates for the source curve 164 * @param srcoff the offset into the array of the beginning of the 165 * the 6 source coordinates 166 * @param left the array for storing the coordinates for the first 167 * half of the subdivided curve 168 * @param leftoff the offset into the array of the beginning of the 169 * the 6 left coordinates 170 * @param right the array for storing the coordinates for the second 171 * half of the subdivided curve 172 * @param rightoff the offset into the array of the beginning of the 173 * the 6 right coordinates 174 * @since 1.2 175 */ 176 public static void subdivide(double src[], int srcoff, 177 double left[], int leftoff, 178 double right[], int rightoff) { 179 double x1 = src[srcoff + 0]; 180 double y1 = src[srcoff + 1]; 181 double ctrlx1 = src[srcoff + 2]; 182 double ctrly1 = src[srcoff + 3]; 183 double ctrlx2 = src[srcoff + 4]; 184 double ctrly2 = src[srcoff + 5]; 185 double x2 = src[srcoff + 6]; 186 double y2 = src[srcoff + 7]; 187 if (left != null) { 188 left[leftoff + 0] = x1; 189 left[leftoff + 1] = y1; 190 } 191 if (right != null) { 192 right[rightoff + 6] = x2; 193 right[rightoff + 7] = y2; 194 } 195 x1 = (x1 + ctrlx1) / 2.0; 196 y1 = (y1 + ctrly1) / 2.0; 197 x2 = (x2 + ctrlx2) / 2.0; 198 y2 = (y2 + ctrly2) / 2.0; 199 double centerx = (ctrlx1 + ctrlx2) / 2.0; 200 double centery = (ctrly1 + ctrly2) / 2.0; 201 ctrlx1 = (x1 + centerx) / 2.0; 202 ctrly1 = (y1 + centery) / 2.0; 203 ctrlx2 = (x2 + centerx) / 2.0; 204 ctrly2 = (y2 + centery) / 2.0; 205 centerx = (ctrlx1 + ctrlx2) / 2.0; 206 centery = (ctrly1 + ctrly2) / 2.0; 207 if (left != null) { 208 left[leftoff + 2] = x1; 209 left[leftoff + 3] = y1; 210 left[leftoff + 4] = ctrlx1; 211 left[leftoff + 5] = ctrly1; 212 left[leftoff + 6] = centerx; 213 left[leftoff + 7] = centery; 214 } 215 if (right != null) { 216 right[rightoff + 0] = centerx; 217 right[rightoff + 1] = centery; 218 right[rightoff + 2] = ctrlx2; 219 right[rightoff + 3] = ctrly2; 220 right[rightoff + 4] = x2; 221 right[rightoff + 5] = y2; 222 } 223 } 224 225 /** 226 * Solves the cubic whose coefficients are in the <code>eqn</code> 227 * array and places the non-complex roots back into the same array, 228 * returning the number of roots. The solved cubic is represented 229 * by the equation: 230 * <pre> 231 * eqn = {c, b, a, d} 232 * dx^3 + ax^2 + bx + c = 0 233 * </pre> 234 * A return value of -1 is used to distinguish a constant equation 235 * that might be always 0 or never 0 from an equation that has no 236 * zeroes. 237 * @param eqn an array containing coefficients for a cubic 238 * @return the number of roots, or -1 if the equation is a constant. 239 * @since 1.2 240 */ 241 public static int solveCubic(double eqn[]) { 242 return solveCubic2(eqn, eqn); 243 } 244 245 /** 246 * Solve the cubic whose coefficients are in the <code>eqn</code> 247 * array and place the non-complex roots into the <code>res</code> 248 * array, returning the number of roots. 249 * The cubic solved is represented by the equation: 250 * eqn = {c, b, a, d} 251 * dx^3 + ax^2 + bx + c = 0 252 * A return value of -1 is used to distinguish a constant equation, 253 * which may be always 0 or never 0, from an equation which has no 254 * zeroes. 255 * @param eqn the specified array of coefficients to use to solve 256 * the cubic equation 257 * @param res the array that contains the non-complex roots 258 * resulting from the solution of the cubic equation 259 * @return the number of roots, or -1 if the equation is a constant 260 * @since 1.3 261 */ 262 public static int solveCubic2(double eqn[], double res[]) { 263 // From Numerical Recipes, 5.6, Quadratic and Cubic Equations 264 double d = eqn[3]; 265 if (d == 0.0) { 266 // The cubic has degenerated to quadratic (or line or ...). 267 return QuadCurve2D.solveQuadratic2(eqn, res); 268 } 269 double a = eqn[2] / d; 270 double b = eqn[1] / d; 271 double c = eqn[0] / d; 272 int roots = 0; 273 double Q = (a * a - 3.0 * b) / 9.0; 274 double R = (2.0 * a * a * a - 9.0 * a * b + 27.0 * c) / 54.0; 275 double R2 = R * R; 276 double Q3 = Q * Q * Q; 277 a = a / 3.0; 278 if (R2 < Q3) { 279 double theta = Math.acos(R / Math.sqrt(Q3)); 280 Q = -2.0 * Math.sqrt(Q); 281 if (res == eqn) { 282 // Copy the eqn so that we don't clobber it with the 283 // roots. This is needed so that fixRoots can do its 284 // work with the original equation. 285 eqn = new double[4]; 286 System.arraycopy(res, 0, eqn, 0, 4); 287 } 288 res[roots++] = Q * Math.cos(theta / 3.0) - a; 289 res[roots++] = Q * Math.cos((theta + Math.PI * 2.0)/ 3.0) - a; 290 res[roots++] = Q * Math.cos((theta - Math.PI * 2.0)/ 3.0) - a; 291 fixRoots(res, eqn); 292 } else { 293 boolean neg = (R < 0.0); 294 double S = Math.sqrt(R2 - Q3); 295 if (neg) { 296 R = -R; 297 } 298 double A = Math.pow(R + S, 1.0 / 3.0); 299 if (!neg) { 300 A = -A; 301 } 302 double B = (A == 0.0) ? 0.0 : (Q / A); 303 res[roots++] = (A + B) - a; 304 } 305 return roots; 306 } 307 308 /* 309 * This pruning step is necessary since solveCubic uses the 310 * cosine function to calculate the roots when there are 3 311 * of them. Since the cosine method can have an error of 312 * +/- 1E-14 we need to make sure that we don't make any 313 * bad decisions due to an error. 314 * 315 * If the root is not near one of the endpoints, then we will 316 * only have a slight inaccuracy in calculating the x intercept 317 * which will only cause a slightly wrong answer for some 318 * points very close to the curve. While the results in that 319 * case are not as accurate as they could be, they are not 320 * disastrously inaccurate either. 321 * 322 * On the other hand, if the error happens near one end of 323 * the curve, then our processing to reject values outside 324 * of the t=[0,1] range will fail and the results of that 325 * failure will be disastrous since for an entire horizontal 326 * range of test points, we will either overcount or undercount 327 * the crossings and get a wrong answer for all of them, even 328 * when they are clearly and obviously inside or outside the 329 * curve. 330 * 331 * To work around this problem, we try a couple of Newton-Raphson 332 * iterations to see if the true root is closer to the endpoint 333 * or further away. If it is further away, then we can stop 334 * since we know we are on the right side of the endpoint. If 335 * we change direction, then either we are now being dragged away 336 * from the endpoint in which case the first condition will cause 337 * us to stop, or we have passed the endpoint and are headed back. 338 * In the second case, we simply evaluate the slope at the 339 * endpoint itself and place ourselves on the appropriate side 340 * of it or on it depending on that result. 341 */ 342 private static void fixRoots(double res[], double eqn[]) { 343 final double EPSILON = 1E-5; 344 for (int i = 0; i < 3; i++) { 345 double t = res[i]; 346 if (Math.abs(t) < EPSILON) { 347 res[i] = findZero(t, 0, eqn); 348 } else if (Math.abs(t - 1) < EPSILON) { 349 res[i] = findZero(t, 1, eqn); 350 } 351 } 352 } 353 354 private static double solveEqn(double eqn[], int order, double t) { 355 double v = eqn[order]; 356 while (--order >= 0) { 357 v = v * t + eqn[order]; 358 } 359 return v; 360 } 361 362 private static double findZero(double t, double target, double eqn[]) { 363 double slopeqn[] = {eqn[1], 2*eqn[2], 3*eqn[3]}; 364 double slope; 365 double origdelta = 0; 366 double origt = t; 367 while (true) { 368 slope = solveEqn(slopeqn, 2, t); 369 if (slope == 0) { 370 // At a local minima - must return 371 return t; 372 } 373 double y = solveEqn(eqn, 3, t); 374 if (y == 0) { 375 // Found it! - return it 376 return t; 377 } 378 // assert(slope != 0 && y != 0); 379 double delta = - (y / slope); 380 // assert(delta != 0); 381 if (origdelta == 0) { 382 origdelta = delta; 383 } 384 if (t < target) { 385 if (delta < 0) return t; 386 } else if (t > target) { 387 if (delta > 0) return t; 388 } else { /* t == target */ 389 return (delta > 0 390 ? (target + java.lang.Double.MIN_VALUE) 391 : (target - java.lang.Double.MIN_VALUE)); 392 } 393 double newt = t + delta; 394 if (t == newt) { 395 // The deltas are so small that we aren't moving... 396 return t; 397 } 398 if (delta * origdelta < 0) { 399 // We have reversed our path. 400 int tag = (origt < t 401 ? getTag(target, origt, t) 402 : getTag(target, t, origt)); 403 if (tag != INSIDE) { 404 // Local minima found away from target - return the middle 405 return (origt + t) / 2; 406 } 407 // Local minima somewhere near target - move to target 408 // and let the slope determine the resulting t. 409 t = target; 410 } else { 411 t = newt; 412 } 413 } 414 } 415 416 /* 417 * Fill an array with the coefficients of the parametric equation 418 * in t, ready for solving against val with solveCubic. 419 * We currently have: 420 * <pre> 421 * val = P(t) = C1(1-t)^3 + 3CP1 t(1-t)^2 + 3CP2 t^2(1-t) + C2 t^3 422 * = C1 - 3C1t + 3C1t^2 - C1t^3 + 423 * 3CP1t - 6CP1t^2 + 3CP1t^3 + 424 * 3CP2t^2 - 3CP2t^3 + 425 * C2t^3 426 * 0 = (C1 - val) + 427 * (3CP1 - 3C1) t + 428 * (3C1 - 6CP1 + 3CP2) t^2 + 429 * (C2 - 3CP2 + 3CP1 - C1) t^3 430 * 0 = C + Bt + At^2 + Dt^3 431 * C = C1 - val 432 * B = 3*CP1 - 3*C1 433 * A = 3*CP2 - 6*CP1 + 3*C1 434 * D = C2 - 3*CP2 + 3*CP1 - C1 435 * </pre> 436 */ 437 private static void fillEqn(double eqn[], double val, 438 double c1, double cp1, double cp2, double c2) { 439 eqn[0] = c1 - val; 440 eqn[1] = (cp1 - c1) * 3.0; 441 eqn[2] = (cp2 - cp1 - cp1 + c1) * 3.0; 442 eqn[3] = c2 + (cp1 - cp2) * 3.0 - c1; 443 } 444 445 private static final int BELOW = -2; 446 private static final int LOWEDGE = -1; 447 private static final int INSIDE = 0; 448 private static final int HIGHEDGE = 1; 449 private static final int ABOVE = 2; 450 451 /* 452 * Determine where coord lies with respect to the range from 453 * low to high. It is assumed that low <= high. The return 454 * value is one of the 5 values BELOW, LOWEDGE, INSIDE, HIGHEDGE, 455 * or ABOVE. 456 */ 457 private static int getTag(double coord, double low, double high) { 458 if (coord <= low) { 459 return (coord < low ? BELOW : LOWEDGE); 460 } 461 if (coord >= high) { 462 return (coord > high ? ABOVE : HIGHEDGE); 463 } 464 return INSIDE; 465 } 466 467 /* 468 * Determine if the pttag represents a coordinate that is already 469 * in its test range, or is on the border with either of the two 470 * opttags representing another coordinate that is "towards the 471 * inside" of that test range. In other words, are either of the 472 * two "opt" points "drawing the pt inward"? 473 */ 474 private static boolean inwards(int pttag, int opt1tag, int opt2tag) { 475 switch (pttag) { 476 case BELOW: 477 case ABOVE: 478 default: 479 return false; 480 case LOWEDGE: 481 return (opt1tag >= INSIDE || opt2tag >= INSIDE); 482 case INSIDE: 483 return true; 484 case HIGHEDGE: 485 return (opt1tag <= INSIDE || opt2tag <= INSIDE); 486 } 487 } 488 489 @Override 490 /** 491 * Creates a new object of the same class as this object. 492 * 493 * @return a clone of this instance. 494 * @exception OutOfMemoryError if there is not enough memory. 495 * @see java.lang.Cloneable 496 * @since 1.2 497 */ 498 public Object clone() { 499 try { 500 return super.clone(); 501 } catch (CloneNotSupportedException e) { 502 // this shouldn't happen, since we are Cloneable 503 throw new InternalError(); 504 } 505 } 506}