cubiccurve2d.java

kaffe Java 解释器语言,源码,Java的子集系统,开放源代码

   * + 6]</code>, its <i>y</i> coordinate at <code>coords[offset +
   * 7]</code>.
   *
   * @param offset the offset of the first coordinate value in
   * <code>coords</code>.
   */
  public static double getFlatness(double[] coords, int offset)
  {
    return Math.sqrt(getFlatnessSq(coords[offset++], coords[offset++],
                                   coords[offset++], coords[offset++],
                                   coords[offset++], coords[offset++],
                                   coords[offset++], coords[offset++]));
  }


  /**
   * Calculates the squared flatness of this curve.  The flatness is
   * the maximal distance of a control point to the line between start
   * and end point.
   *
   * <p><img src="doc-files/CubicCurve2D-4.png" width="350" height="180"
   * alt="A drawing that illustrates the flatness" />
   *
   * <p>In the above drawing, the straight line connecting start point
   * P1 and end point P2 is depicted in gray.  In comparison to C1,
   * control point C2 is father away from the gray line. Therefore,
   * the result will be the square of the distance between C2 and the
   * gray line, i.e. the squared length of the red line.
   */
  public double getFlatnessSq()
  {
    return getFlatnessSq(getX1(), getY1(), getCtrlX1(), getCtrlY1(),
                         getCtrlX2(), getCtrlY2(), getX2(), getY2());
  }


  /**
   * Calculates the flatness of this curve.  The flatness is the
   * maximal distance of a control point to the line between start and
   * end point.
   *
   * <p><img src="doc-files/CubicCurve2D-4.png" width="350" height="180"
   * alt="A drawing that illustrates the flatness" />
   *
   * <p>In the above drawing, the straight line connecting start point
   * P1 and end point P2 is depicted in gray.  In comparison to C1,
   * control point C2 is father away from the gray line. Therefore,
   * the result will be the distance between C2 and the gray line,
   * i.e. the length of the red line.
   */
  public double getFlatness()
  {
    return Math.sqrt(getFlatnessSq(getX1(), getY1(), getCtrlX1(),
                                   getCtrlY1(), getCtrlX2(), getCtrlY2(),
                                   getX2(), getY2()));
  }


  /**
   * Subdivides this curve into two halves.
   *
   * <p><img src="doc-files/CubicCurve2D-3.png" width="700"
   * height="180" alt="A drawing that illustrates the effects of
   * subdividing a CubicCurve2D" />
   *
   * @param left a curve whose geometry will be set to the left half
   * of this curve, or <code>null</code> if the caller is not
   * interested in the left half.
   *
   * @param right a curve whose geometry will be set to the right half
   * of this curve, or <code>null</code> if the caller is not
   * interested in the right half.
   */
  public void subdivide(CubicCurve2D left, CubicCurve2D right)
  {
    // Use empty slots at end to share single array.
    double[] d = new double[] { getX1(), getY1(), getCtrlX1(), getCtrlY1(),
                                getCtrlX2(), getCtrlY2(), getX2(), getY2(),
                                0, 0, 0, 0, 0, 0 };
    subdivide(d, 0, d, 0, d, 6);
    if (left != null)
      left.setCurve(d, 0);
    if (right != null)
      right.setCurve(d, 6);
  }


  /**
   * Subdivides a cubic curve into two halves.
   *
   * <p><img src="doc-files/CubicCurve2D-3.png" width="700"
   * height="180" alt="A drawing that illustrates the effects of
   * subdividing a CubicCurve2D" />
   *
   * @param src the curve to be subdivided.
   *
   * @param left a curve whose geometry will be set to the left half
   * of <code>src</code>, or <code>null</code> if the caller is not
   * interested in the left half.
   *
   * @param right a curve whose geometry will be set to the right half
   * of <code>src</code>, or <code>null</code> if the caller is not
   * interested in the right half.
   */
  public static void subdivide(CubicCurve2D src,
                               CubicCurve2D left, CubicCurve2D right)
  {
    src.subdivide(left, right);
  }


  /**
   * Subdivides a cubic curve into two halves, passing all coordinates
   * in an array.
   *
   * <p><img src="doc-files/CubicCurve2D-3.png" width="700"
   * height="180" alt="A drawing that illustrates the effects of
   * subdividing a CubicCurve2D" />
   *
   * <p>The left end point and the right start point will always be
   * identical. Memory-concious programmers thus may want to pass the
   * same array for both <code>left</code> and <code>right</code>, and
   * set <code>rightOff</code> to <code>leftOff + 6</code>.
   *
   * @param src an array containing the coordinates of the curve to be
   * subdivided.  The <i>x</i> coordinate of the start point P1 is
   * located at <code>src[srcOff]</code>, its <i>y</i> at
   * <code>src[srcOff + 1]</code>.  The <i>x</i> coordinate of the
   * first control point C1 is located at <code>src[srcOff +
   * 2]</code>, its <i>y</i> at <code>src[srcOff + 3]</code>.  The
   * <i>x</i> coordinate of the second control point C2 is located at
   * <code>src[srcOff + 4]</code>, its <i>y</i> at <code>src[srcOff +
   * 5]</code>. The <i>x</i> coordinate of the end point is located at
   * <code>src[srcOff + 6]</code>, its <i>y</i> at <code>src[srcOff +
   * 7]</code>.
   *
   * @param srcOff an offset into <code>src</code>, specifying
   * the index of the start point&#x2019;s <i>x</i> coordinate.
   *
   * @param left an array that will receive the coordinates of the
   * left half of <code>src</code>. It is acceptable to pass
   * <code>src</code>. A caller who is not interested in the left half
   * can pass <code>null</code>.
   *
   * @param leftOff an offset into <code>left</code>, specifying the
   * index where the start point&#x2019;s <i>x</i> coordinate will be
   * stored.
   *
   * @param right an array that will receive the coordinates of the
   * right half of <code>src</code>. It is acceptable to pass
   * <code>src</code> or <code>left</code>. A caller who is not
   * interested in the right half can pass <code>null</code>.
   *
   * @param rightOff an offset into <code>right</code>, specifying the
   * index where the start point&#x2019;s <i>x</i> coordinate will be
   * stored.
   */
  public static void subdivide(double[] src, int srcOff,
                               double[] left, int leftOff,
                               double[] right, int rightOff)
  {
    // To understand this code, please have a look at the image
    // "CubicCurve2D-3.png" in the sub-directory "doc-files".
    double src_C1_x, src_C1_y, src_C2_x, src_C2_y;
    double left_P1_x, left_P1_y;
    double left_C1_x, left_C1_y, left_C2_x, left_C2_y;
    double right_C1_x, right_C1_y, right_C2_x, right_C2_y;
    double right_P2_x, right_P2_y;
    double Mid_x, Mid_y; // Mid = left.P2 = right.P1

    left_P1_x = src[srcOff];
    left_P1_y = src[srcOff + 1];
    src_C1_x = src[srcOff + 2];
    src_C1_y = src[srcOff + 3];
    src_C2_x = src[srcOff + 4];
    src_C2_y = src[srcOff + 5];
    right_P2_x = src[srcOff + 6];
    right_P2_y = src[srcOff + 7];

    left_C1_x = (left_P1_x + src_C1_x) / 2;
    left_C1_y = (left_P1_y + src_C1_y) / 2;
    right_C2_x = (right_P2_x + src_C2_x) / 2;
    right_C2_y = (right_P2_y + src_C2_y) / 2;
    Mid_x = (src_C1_x + src_C2_x) / 2;
    Mid_y = (src_C1_y + src_C2_y) / 2;
    left_C2_x = (left_C1_x + Mid_x) / 2;
    left_C2_y = (left_C1_y + Mid_y) / 2;
    right_C1_x = (Mid_x + right_C2_x) / 2;
    right_C1_y = (Mid_y + right_C2_y) / 2;
    Mid_x = (left_C2_x + right_C1_x) / 2;
    Mid_y = (left_C2_y + right_C1_y) / 2;

    if (left != null)
    {
      left[leftOff] = left_P1_x;
      left[leftOff + 1] = left_P1_y;
      left[leftOff + 2] = left_C1_x;
      left[leftOff + 3] = left_C1_y;
      left[leftOff + 4] = left_C2_x;
      left[leftOff + 5] = left_C2_y;
      left[leftOff + 6] = Mid_x;
      left[leftOff + 7] = Mid_y;
    }

    if (right != null)
    {
      right[rightOff] = Mid_x;
      right[rightOff + 1] = Mid_y;
      right[rightOff + 2] = right_C1_x;
      right[rightOff + 3] = right_C1_y;
      right[rightOff + 4] = right_C2_x;
      right[rightOff + 5] = right_C2_y;
      right[rightOff + 6] = right_P2_x;
      right[rightOff + 7] = right_P2_y;
    }
  }


  /**
   * Finds the non-complex roots of a cubic equation, placing the
   * results into the same array as the equation coefficients. The
   * following equation is being solved:
   *
   * <blockquote><code>eqn[3]</code> &#xb7; <i>x</i><sup>3</sup>
   * + <code>eqn[2]</code> &#xb7; <i>x</i><sup>2</sup>
   * + <code>eqn[1]</code> &#xb7; <i>x</i>
   * + <code>eqn[0]</code>
   * = 0
   * </blockquote>
   *
   * <p>For some background about solving cubic equations, see the
   * article <a
   * href="http://planetmath.org/encyclopedia/CubicFormula.html"
   * >&#x201c;Cubic Formula&#x201d;</a> in <a
   * href="http://planetmath.org/" >PlanetMath</a>.  For an extensive
   * library of numerical algorithms written in the C programming
   * language, see the <a href= "http://www.gnu.org/software/gsl/">GNU
   * Scientific Library</a>, from which this implementation was
   * adapted.
   *
   * @param eqn an array with the coefficients of the equation. When
   * this procedure has returned, <code>eqn</code> will contain the
   * non-complex solutions of the equation, in no particular order.
   *
   * @return the number of non-complex solutions. A result of 0
   * indicates that the equation has no non-complex solutions. A
   * result of -1 indicates that the equation is constant (i.e.,
   * always or never zero).
   *
   * @see #solveCubic(double[], double[])
   * @see QuadCurve2D#solveQuadratic(double[],double[])
   *
   * @author <a href="mailto:[e-mail]">Brian Gough</a>
   * (original C implementation in the <a href=
   * "http://www.gnu.org/software/gsl/">GNU Scientific Library</a>)
   *
   * @author <a href="mailto:[e-mail]">Sascha Brawer</a>
   * (adaptation to Java)
   */
  public static int solveCubic(double[] eqn)
  {
    return solveCubic(eqn, eqn);
  }


  /**
   * Finds the non-complex roots of a cubic equation. The following
   * equation is being solved:
   *
   * <blockquote><code>eqn[3]</code> &#xb7; <i>x</i><sup>3</sup>
   * + <code>eqn[2]</code> &#xb7; <i>x</i><sup>2</sup>
   * + <code>eqn[1]</code> &#xb7; <i>x</i>
   * + <code>eqn[0]</code>
   * = 0
   * </blockquote>
   *
   * <p>For some background about solving cubic equations, see the
   * article <a
   * href="http://planetmath.org/encyclopedia/CubicFormula.html"
   * >&#x201c;Cubic Formula&#x201d;</a> in <a
   * href="http://planetmath.org/" >PlanetMath</a>.  For an extensive
   * library of numerical algorithms written in the C programming
   * language, see the <a href= "http://www.gnu.org/software/gsl/">GNU
   * Scientific Library</a>, from which this implementation was
   * adapted.
   *
   * @see QuadCurve2D#solveQuadratic(double[],double[])
   *
   * @param eqn an array with the coefficients of the equation.
   *
   * @param res an array into which the non-complex roots will be
   * stored.  The results may be in an arbitrary order. It is safe to
   * pass the same array object reference for both <code>eqn</code>
   * and <code>res</code>.
   *
   * @return the number of non-complex solutions. A result of 0
   * indicates that the equation has no non-complex solutions. A
   * result of -1 indicates that the equation is constant (i.e.,
   * always or never zero).
   *
   * @author <a href="mailto:[e-mail]">Brian Gough</a>
   * (original C implementation in the <a href=
   * "http://www.gnu.org/software/gsl/">GNU Scientific Library</a>)
   *
   * @author <a href="mailto:[e-mail]">Sascha Brawer</a>
   * (adaptation to Java)
   */
  public static int solveCubic(double[] eqn, double[] res)
  {
    // Adapted from poly/solve_cubic.c in the GNU Scientific Library
    // (GSL), revision 1.7 of 2003-07-26. For the original source, see
    // http://www.gnu.org/software/gsl/
    //
    // Brian Gough, the author of that code, has granted the
    // permission to use it in GNU Classpath under the GNU Classpath
    // license, and has assigned the copyright to the Free Software
    // Foundation.
    //
    // The Java implementation is very similar to the GSL code, but
    // not a strict one-to-one copy. For example, GSL would sort the
    // result.

    double a, b, c, q, r, Q, R;
    double c3, Q3, R2, CR2, CQ3;

    // If the cubic coefficient is zero, we have a quadratic equation.
    c3 = eqn[3];
    if (c3 == 0)
      return QuadCurve2D.solveQuadratic(eqn, res);

    // Divide the equation by the cubic coefficient.
    c = eqn[0] / c3;
    b = eqn[1] / c3;
    a = eqn[2] / c3;

    // We now need to solve x^3 + ax^2 + bx + c = 0.
    q = a * a - 3 * b;
    r = 2 * a * a * a - 9 * a * b + 27 * c;

    Q = q / 9;
    R = r / 54;

    Q3 = Q * Q * Q;
    R2 = R * R;

    CR2 = 729 * r * r;
    CQ3 = 2916 * q * q * q;

    if (R == 0 && Q == 0)
    {
      // The GNU Scientific Library would return three identical
      // solutions in this case.
      res[0] = -a/3;
      return 1;
    }

    if (CR2 == CQ3) 
    {
      /* this test is actually R2 == Q3, written in a form suitable
         for exact computation with integers */

      /* Due to finite precision some double roots may be missed, and
         considered to be a pair of complex roots z = x +/- epsilon i
         close to the real axis. */

      double sqrtQ = Math.sqrt(Q);

      if (R > 0)
      {
        res[0] = -2 * sqrtQ - a/3;
        res[1] = sqrtQ - a/3;
      }
      else
      {
        res[0] = -sqrtQ - a/3;
        res[1] = 2 * sqrtQ - a/3;
      }
      return 2;
    }

    if (CR2 < CQ3) /* equivalent to R2 < Q3 */
    {
      double sqrtQ = Math.sqrt(Q);
      double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
      double theta = Math.acos(R / sqrtQ3);
      double norm = -2 * sqrtQ;
      res[0] = norm * Math.cos(theta / 3) - a / 3;
      res[1] = norm * Math.cos((theta + 2.0 * Math.PI) / 3) - a/3;
      res[2] = norm * Math.cos((theta - 2.0 * Math.PI) / 3) - a/3;

      // The GNU Scientific Library sorts the results. We don't.
      return 3;
    }

    double sgnR = (R >= 0 ? 1 : -1);
    double A = -sgnR * Math.pow(Math.abs(R) + Math.sqrt(R2 - Q3), 1.0/3.0);
    double B = Q / A ;
    res[0] = A + B - a/3;
    return 1;
  }


  /**
   * Determines whether a position lies inside the area that is bounded
   * by the curve and the straight line connecting its end points.
   *
   * <p><img src="doc-files/CubicCurve2D-5.png" width="350" height="180"
   * alt="A drawing of the area spanned by the curve" />