Bezier.java

  1. package org.opentrafficsim.core.geometry;

  2. import java.awt.geom.Line2D;

  3. import org.djutils.draw.point.OrientedPoint2d;
  4. import org.djutils.draw.point.Point2d;
  5. import org.djutils.exceptions.Throw;
  6. import org.opentrafficsim.base.geometry.OtsGeometryException;
  7. import org.opentrafficsim.base.geometry.OtsLine2d;

  8. /**
  9.  * Generation of B&eacute;zier curves. <br>
  10.  * The class implements the cubic(...) method to generate a cubic B&eacute;zier curve using the following formula: B(t) = (1 -
  11.  * t)<sup>3</sup>P<sub>0</sub> + 3t(1 - t)<sup>2</sup> P<sub>1</sub> + 3t<sup>2</sup> (1 - t) P<sub>2</sub> + t<sup>3</sup>
  12.  * P<sub>3</sub> where P<sub>0</sub> and P<sub>3</sub> are the end points, and P<sub>1</sub> and P<sub>2</sub> the control
  13.  * points. <br>
  14.  * For a smooth movement, one of the standard implementations of the cubic(...) function offered is the case where P<sub>1</sub>
  15.  * is positioned halfway between P<sub>0</sub> and P<sub>3</sub> starting from P<sub>0</sub> in the direction of P<sub>3</sub>,
  16.  * and P<sub>2</sub> is positioned halfway between P<sub>3</sub> and P<sub>0</sub> starting from P<sub>3</sub> in the direction
  17.  * of P<sub>0</sub>.<br>
  18.  * Finally, an n-point generalization of the B&eacute;zier curve is implemented with the bezier(...) function.
  19.  * <p>
  20.  * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
  21.  * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
  22.  * </p>
  23.  * @author <a href="https://github.com/averbraeck">Alexander Verbraeck</a>
  24.  * @author <a href="https://github.com/peter-knoppers">Peter Knoppers</a>
  25.  */
  26. public final class Bezier
  27. {
  28.     /** The default number of points to use to construct a B&eacute;zier curve. */
  29.     private static final int DEFAULT_NUM_POINTS = 64;

  30.     /** Cached factorial values. */
  31.     private static long[] fact = new long[] {1L, 1L, 2L, 6L, 24L, 120L, 720L, 5040L, 40320L, 362880L, 3628800L, 39916800L,
  32.             479001600L, 6227020800L, 87178291200L, 1307674368000L, 20922789888000L, 355687428096000L, 6402373705728000L,
  33.             121645100408832000L, 2432902008176640000L};

  34.     /** Utility class. */
  35.     private Bezier()
  36.     {
  37.         // do not instantiate
  38.     }

  39.     /**
  40.      * Construct a cubic B&eacute;zier curve from start to end with two control points.
  41.      * @param numPoints the number of points for the B&eacute;zier curve
  42.      * @param start the start point of the B&eacute;zier curve
  43.      * @param control1 the first control point
  44.      * @param control2 the second control point
  45.      * @param end the end point of the B&eacute;zier curve
  46.      * @return a cubic B&eacute;zier curve between start and end, with the two provided control points
  47.      * @throws OtsGeometryException in case the number of points is less than 2 or the B&eacute;zier curve could not be
  48.      *             constructed
  49.      */
  50.     public static OtsLine2d cubic(final int numPoints, final Point2d start, final Point2d control1, final Point2d control2,
  51.             final Point2d end) throws OtsGeometryException
  52.     {
  53.         Throw.when(numPoints < 2, OtsGeometryException.class, "Number of points too small (got %d; minimum value is 2)",
  54.                 numPoints);
  55.         Point2d[] points = new Point2d[numPoints];
  56.         for (int n = 0; n < numPoints; n++)
  57.         {
  58.             double t = n / (numPoints - 1.0);
  59.             double x = B3(t, start.x, control1.x, control2.x, end.x);
  60.             double y = B3(t, start.y, control1.y, control2.y, end.y);
  61.             points[n] = new Point2d(x, y);
  62.         }
  63.         return new OtsLine2d(points);
  64.     }

  65.     /**
  66.      * Construct a cubic B&eacute;zier curve from start to end with two generated control points at half the distance between
  67.      * start and end. The z-value is interpolated in a linear way.
  68.      * @param numPoints the number of points for the B&eacute;zier curve
  69.      * @param start the directed start point of the B&eacute;zier curve
  70.      * @param end the directed end point of the B&eacute;zier curve
  71.      * @return a cubic B&eacute;zier curve between start and end, with the two provided control points
  72.      * @throws OtsGeometryException in case the number of points is less than 2 or the B&eacute;zier curve could not be
  73.      *             constructed
  74.      */
  75.     public static OtsLine2d cubic(final int numPoints, final OrientedPoint2d start, final OrientedPoint2d end)
  76.             throws OtsGeometryException
  77.     {
  78.         return cubic(numPoints, start, end, 1.0);
  79.     }

  80.     /**
  81.      * Construct a cubic B&eacute;zier curve from start to end with two generated control points at half the distance between
  82.      * start and end. The z-value is interpolated in a linear way.
  83.      * @param numPoints the number of points for the B&eacute;zier curve
  84.      * @param start the directed start point of the B&eacute;zier curve
  85.      * @param end the directed end point of the B&eacute;zier curve
  86.      * @param shape 1 = control points at half the distance between start and end, &gt; 1 results in a pointier shape, &lt; 1
  87.      *            results in a flatter shape, value should be above 0
  88.      * @return a cubic B&eacute;zier curve between start and end, with the two determined control points
  89.      * @throws OtsGeometryException in case the number of points is less than 2 or the B&eacute;zier curve could not be
  90.      *             constructed
  91.      */
  92.     public static OtsLine2d cubic(final int numPoints, final OrientedPoint2d start, final OrientedPoint2d end,
  93.             final double shape) throws OtsGeometryException
  94.     {
  95.         return cubic(numPoints, start, end, shape, false);
  96.     }

  97.     /**
  98.      * Construct a cubic B&eacute;zier curve from start to end with two generated control points at half the distance between
  99.      * start and end. The z-value is interpolated in a linear way.
  100.      * @param numPoints the number of points for the B&eacute;zier curve
  101.      * @param start the directed start point of the B&eacute;zier curve
  102.      * @param end the directed end point of the B&eacute;zier curve
  103.      * @param shape shape factor; 1 = control points at half the distance between start and end, &gt; 1 results in a pointier
  104.      *            shape, &lt; 1 results in a flatter shape, value should be above 0
  105.      * @param weighted control point distance relates to distance to projected point on extended line from other end
  106.      * @return a cubic B&eacute;zier curve between start and end, with the two determined control points
  107.      * @throws OtsGeometryException in case the number of points is less than 2 or the B&eacute;zier curve could not be
  108.      *             constructed
  109.      */
  110.     public static OtsLine2d cubic(final int numPoints, final OrientedPoint2d start, final OrientedPoint2d end,
  111.             final double shape, final boolean weighted) throws OtsGeometryException
  112.     {
  113.         return bezier(cubicControlPoints(start, end, shape, weighted));
  114.     }

  115.     /**
  116.      * Construct control points for a cubic B&eacute;zier curve from start to end with two generated control points at half the
  117.      * distance between start and end.
  118.      * @param start the directed start point of the B&eacute;zier curve
  119.      * @param end the directed end point of the B&eacute;zier curve
  120.      * @param shape shape factor; 1 = control points at half the distance between start and end, &gt; 1 results in a pointier
  121.      *            shape, &lt; 1 results in a flatter shape, value should be above 0
  122.      * @param weighted control point distance relates to distance to projected point on extended line from other end
  123.      * @return a cubic B&eacute;zier curve between start and end, with the two determined control points
  124.      */
  125.     public static Point2d[] cubicControlPoints(final OrientedPoint2d start, final OrientedPoint2d end, final double shape,
  126.             final boolean weighted)
  127.     {
  128.         Throw.when(shape <= 0.0, IllegalArgumentException.class, "Shape factor must be above 0.0.");
  129.         Point2d control1;
  130.         Point2d control2;

  131.         if (weighted)
  132.         {
  133.             // each control point is 'w' * the distance between the end-points away from the respective end point
  134.             // 'w' is a weight given by the distance from the end point to the extended line of the other end point
  135.             double dx = end.x - start.x;
  136.             double dy = end.y - start.y;
  137.             double distance = shape * Math.hypot(dx, dy);
  138.             double cosEnd = Math.cos(end.getDirZ());
  139.             double sinEnd = Math.sin(end.getDirZ());
  140.             double dStart = Line2D.ptLineDist(end.x, end.y, end.x + cosEnd, end.y + sinEnd, start.x, start.y);
  141.             double cosStart = Math.cos(start.getDirZ());
  142.             double sinStart = Math.sin(start.getDirZ());
  143.             double dEnd = Line2D.ptLineDist(start.x, start.y, start.x + cosStart, start.y + sinStart, end.x, end.y);
  144.             double wStart = dStart / (dStart + dEnd);
  145.             double wEnd = dEnd / (dStart + dEnd);
  146.             double wStartDistance = wStart * distance;
  147.             double wEndDistance = wEnd * distance;
  148.             control1 = new Point2d(start.x + wStartDistance * cosStart, start.y + wStartDistance * sinStart);
  149.             // - (minus) as the angle is where the line leaves, i.e. from shape point to end
  150.             control2 = new Point2d(end.x - wEndDistance * cosEnd, end.y - wEndDistance * sinEnd);
  151.         }
  152.         else
  153.         {
  154.             // each control point is half the distance between the end-points away from the respective end point
  155.             double dx = end.x - start.x;
  156.             double dy = end.y - start.y;
  157.             double distance2 = shape * .5 * Math.hypot(dx, dy);
  158.             control1 = new Point2d(start.x + distance2 * Math.cos(start.getDirZ()),
  159.                     start.y + distance2 * Math.sin(start.getDirZ()));
  160.             control2 = new Point2d(end.x - distance2 * Math.cos(end.getDirZ()), end.y - distance2 * Math.sin(end.getDirZ()));
  161.         }

  162.         return new Point2d[] {start, control1, control2, end};
  163.     }

  164.     /**
  165.      * Construct a cubic B&eacute;zier curve from start to end with two generated control points at half the distance between
  166.      * start and end. The z-value is interpolated in a linear way.
  167.      * @param start the directed start point of the B&eacute;zier curve
  168.      * @param end the directed end point of the B&eacute;zier curve
  169.      * @return a cubic B&eacute;zier curve between start and end, with the two provided control points
  170.      * @throws OtsGeometryException in case the number of points is less than 2 or the B&eacute;zier curve could not be
  171.      *             constructed
  172.      */
  173.     public static OtsLine2d cubic(final OrientedPoint2d start, final OrientedPoint2d end) throws OtsGeometryException
  174.     {
  175.         return cubic(DEFAULT_NUM_POINTS, start, end);
  176.     }

  177.     /**
  178.      * Calculate the cubic B&eacute;zier point with B(t) = (1 - t)<sup>3</sup>P<sub>0</sub> + 3t(1 - t)<sup>2</sup>
  179.      * P<sub>1</sub> + 3t<sup>2</sup> (1 - t) P<sub>2</sub> + t<sup>3</sup> P<sub>3</sub>.
  180.      * @param t the fraction
  181.      * @param p0 the first point of the curve
  182.      * @param p1 the first control point
  183.      * @param p2 the second control point
  184.      * @param p3 the end point of the curve
  185.      * @return the cubic bezier value B(t)
  186.      */
  187.     @SuppressWarnings("checkstyle:methodname")
  188.     private static double B3(final double t, final double p0, final double p1, final double p2, final double p3)
  189.     {
  190.         double t2 = t * t;
  191.         double t3 = t2 * t;
  192.         double m = (1.0 - t);
  193.         double m2 = m * m;
  194.         double m3 = m2 * m;
  195.         return m3 * p0 + 3.0 * t * m2 * p1 + 3.0 * t2 * m * p2 + t3 * p3;
  196.     }

  197.     /**
  198.      * Construct a B&eacute;zier curve of degree n.
  199.      * @param numPoints the number of points for the B&eacute;zier curve to be constructed
  200.      * @param points the points of the curve, where the first and last are begin and end point, and the intermediate ones are
  201.      *            control points. There should be at least two points.
  202.      * @return the B&eacute;zier value B(t) of degree n, where n is the number of points in the array
  203.      * @throws OtsGeometryException in case the number of points is less than 2 or the B&eacute;zier curve could not be
  204.      *             constructed
  205.      */
  206.     public static OtsLine2d bezier(final int numPoints, final Point2d... points) throws OtsGeometryException
  207.     {
  208.         Point2d[] result = new Point2d[numPoints];
  209.         double[] px = new double[points.length];
  210.         double[] py = new double[points.length];
  211.         for (int i = 0; i < points.length; i++)
  212.         {
  213.             px[i] = points[i].x;
  214.             py[i] = points[i].y;
  215.         }
  216.         for (int n = 0; n < numPoints; n++)
  217.         {
  218.             double t = n / (numPoints - 1.0);
  219.             double x = Bn(t, px);
  220.             double y = Bn(t, py);
  221.             result[n] = new Point2d(x, y);
  222.         }
  223.         return new OtsLine2d(result);
  224.     }

  225.     /**
  226.      * Construct a B&eacute;zier curve of degree n.
  227.      * @param points the points of the curve, where the first and last are begin and end point, and the intermediate ones are
  228.      *            control points. There should be at least two points.
  229.      * @return the B&eacute;zier value B(t) of degree n, where n is the number of points in the array
  230.      * @throws OtsGeometryException in case the number of points is less than 2 or the B&eacute;zier curve could not be
  231.      *             constructed
  232.      */
  233.     public static OtsLine2d bezier(final Point2d... points) throws OtsGeometryException
  234.     {
  235.         return bezier(DEFAULT_NUM_POINTS, points);
  236.     }

  237.     /**
  238.      * Calculate the B&eacute;zier point of degree n, with B(t) = Sum(i = 0..n) [C(n, i) * (1 - t)<sup>n-i</sup> t<sup>i</sup>
  239.      * P<sub>i</sub>], where C(n, k) is the binomial coefficient defined by n! / ( k! (n-k)! ), ! being the factorial operator.
  240.      * @param t the fraction
  241.      * @param p the points of the curve, where the first and last are begin and end point, and the intermediate ones are control
  242.      *            points
  243.      * @return the B&eacute;zier value B(t) of degree n, where n is the number of points in the array
  244.      */
  245.     @SuppressWarnings("checkstyle:methodname")
  246.     static double Bn(final double t, final double... p)
  247.     {
  248.         double b = 0.0;
  249.         double m = (1.0 - t);
  250.         int n = p.length - 1;
  251.         double fn = factorial(n);
  252.         for (int i = 0; i <= n; i++)
  253.         {
  254.             double c = fn / (factorial(i) * (factorial(n - i)));
  255.             b += c * Math.pow(m, n - i) * Math.pow(t, i) * p[i];
  256.         }
  257.         return b;
  258.     }

  259.     /**
  260.      * Calculate factorial(k), which is k * (k-1) * (k-2) * ... * 1. For factorials up to 20, a lookup table is used.
  261.      * @param k the parameter
  262.      * @return factorial(k)
  263.      */
  264.     private static double factorial(final int k)
  265.     {
  266.         if (k < fact.length)
  267.         {
  268.             return fact[k];
  269.         }
  270.         double f = 1;
  271.         for (int i = 2; i <= k; i++)
  272.         {
  273.             f = f * i;
  274.         }
  275.         return f;
  276.     }

  277. }