To print the points on the Hilbert curve, the printPoints() method in Listing C.18 starts with a point at the origin \((0,0)\) and repeatedly updates its position using the directions computed by hilbertCurve().
Listing C.18ch4 / HilbertCurve
// Print the coordinates of all points along a Hilbert curve.
private static void printPoints(int order) {
var directions = new ArrayList<Integer>();
hilbertCurve(order, NORTH, +1, directions);
Point p = new Point(0, 0);
for (int d : directions) {
System.out.printf("(%d,%d) ", (int) p.x(), (int) p.y());
p = switch (d) {
case NORTH -> new Point(p.x(), p.y() + 1);
case EAST -> new Point(p.x() + 1, p.y());
case SOUTH -> new Point(p.x(), p.y() - 1);
case WEST -> new Point(p.x() - 1, p.y());
default -> throw new RuntimeException("Invalid direction");
};
}
System.out.printf("(%d,%d)\n", (int) p.x(), (int) p.y());
}
The subdivide() method shown in Listing C.20 draws one edge of the Koch snowflake. It takes the end points of a line segment as its input and recursively subdivides it into a fractal curve. We first split the segment between
start and end into three equal parts and then erect an equilateral triangle above the two inner points p and q. The formula for the tip of this triangle tip can be derived using basic trigonometry. The resulting four line
segments are then subdivided further.
Listing C.20ch4 / KochCurve
// Draw one edge of the Koch snowflake.
public static void subdivide(Point start, Point end, int depth,
Graphics g) {
if (depth == 0) {
g.drawLine((int) start.x(), (int) start.y(),
(int) end.x(), (int) end.y());
} else {
Point p = Point.lerp(start, end, 1.0 / 3.0);
Point q = Point.lerp(start, end, 2.0 / 3.0);
Point tip = new Point(
(q.x() + p.x()) / 2 - (q.y() - p.y()) * Math.sqrt(3) / 2,
(q.y() + p.y()) / 2 + (q.x() - p.x()) * Math.sqrt(3) / 2);
subdivide(start, p, depth - 1, g);
subdivide(p, tip, depth - 1, g);
subdivide(tip, q, depth - 1, g);
subdivide(q, end, depth - 1, g);
}
}
b)The implementation of LineSegment is shown in Listing C.21. To implement pointAt(), we use the lerp() function to compute a point between the two end points. To implement closestPoint(), we first determine the parameter of the point closest to m along the curve using closestParam() and then return the point for that parameter. The closestParam() method uses Eq. (C.4) to find the parameter for the closest point along a line and then clamps it to the range \([0,1]\).
Listing C.21ch4 / LineSegment
// A data type for representing line segments in the plane.
public record LineSegment(Point p0, Point p1) {
public Point pointAt(double t) {
return Point.lerp(p0, p1, t);
}
public double closestParam(Point m) {
double d = Point.distance(p0, p1);
if (d == 0.0)
return 0.0;
double t = ((m.x() - p0.x()) * (p1.x() - p0.x())
+ (m.y() - p0.y()) * (p1.y() - p0.y()))
/ (d * d);
return Math.max(0.0f, Math.min(1.0f, t));
}
public Point closestPoint(Point m) {
return pointAt(closestParam(m));
}
}
c)To find the closest point on a Bézier curve, we apply de Casteljau’s algorithm as shown in Listing C.22. If the current segment is sufficiently straight, we estimate the closest point and its distance from m. Otherwise, we split the segment into two curves, recursively find the closest point on each, and return the point that is closer to m.
Listing C.22ch4 / BezierCurve
public record PointAndDist(Point point, double distance) {
}
// Compute the point on the curve that is closest to 'm'.
public PointAndDist closestPoint(Point m, double tolerance) {
if (isStraight(tolerance)) {
double t = new LineSegment(p0(), p3()).closestParam(m);
Point p = pointAt(t);
return new PointAndDist(p, Point.distance(p, m));
} else {
var segments = split();
var left = segments[0].closestPoint(m, tolerance);
var right = segments[1].closestPoint(m, tolerance);
return left.distance < right.distance ? left : right;
}
}
b)It’s not hard to prove that the formula for a quadratic Bézier curve simplifies to a line segment if the control point \(\vec {p}_1\) lies halfway between the end points \(\vec {p}_0\) and \(\vec {p}_2\). A simple straightness test is to compare the distance of \(\vec {p}_1\) to the
midpoint; this is implemented in isStraight() in Listing C.23.
Listing C.23ch4 / QuadraticBezier
public record QuadraticBezier(Point p0, Point p1, Point p2) {
private boolean isStraight(double tolerance) {
return Point.distance(p1, Point.lerp(p0, p2, 0.5)) < tolerance;
}
public void subdivide(double tolerance, List<Point> points) {
if (isStraight(tolerance))
points.add(p2);
Point q0 = Point.lerp(p0, p1, 0.5);
Point q1 = Point.lerp(p1, p2, 0.5);
Point s = Point.lerp(q0, q1, 0.5);
new QuadraticBezier(p0, q0, s).subdivide(tolerance, points);
new QuadraticBezier(s, q1, p2).subdivide(tolerance, points);
}
}
c)As indicated in Fig. 4.10, each step of the construction yields one point on the curve and two new control points, one for each half of the original curve. This give us the subdivision scheme implemented by the subdivide() method in Listing C.23.
Solution to Exercise 4.8
a)From Eq. (4.2) we have \(\vec {B}(1/2) = (\vec {p}_0 + 3\vec {p}_1 + 3\vec {p}_2 + \vec {p}_3)/8\). The sequence of midpoints computed by de Casteljau’s
algorithm computes the same point:
b)If you construct a diagram similar to the one in part (a) but combine adjacent points \(\vec {a}\) and \(\vec {b}\) in each row using linear interpolation \((1-t)\vec {a}+t\vec {b}\), you obtain the formula for Bézier curves Eq. (4.2) on the final line.
c)If we set \(\vec {p}_1=\lerp (\vec {p}_0, \vec {p}_3, 1/3)\) and \(\vec {p}_2=\lerp (\vec {p}_0, \vec {p}_3, 2/3)\) and plug these values into the definition of the Bézier curves, we obtain:
\(\seteqnumber{0}{C.}{4}\)
\begin{align*}
\vec {B}(t) &= \left [(1-t)^3 + 2t(1-t)^2+t^2(1-t)\right ]\vec {p}_0 \\&\qquad + \left [t(1-t)^2 + 2t^2(1-t)+t^3\right ]\vec {p}_3\\ &= (1-t)[1-t+t]^2\vec {p}_0 + t[t+1-t]^2\vec {p}_3\\ &= (1-t)\vec {p}_0 + t\vec {p}_3.
\end{align*}
This is the line segment from \(\vec {p}_0\) to \(\vec {p}_3\).
d)According to Eq. (4.2), every point \(\vec {B}(t)\) on the Bézier curve is a linear combination of the four control points
In the second to last step we have used the fact that the weights \(w_k\) sum to 1 (Eq. (4.3)).
Solution to Exercise 4.9
a)The five-pointed star can be constructed in two ways, either by joining the sides with single arrows or the sides with double arrows. In the former case the arrows point outwards, in the latter inwards. Placing even two thin rhombs side by side isn’t allowed because there is no way to
properly line up the arrow symbols, so there is no way to form a ten-pointed star.
b)Deflating a thick rhomb twice produces three thick rhombs, three thin rhombs, and four halves of a thick rhomb:
c)Listing C.24 defines two functions deflateThin() and deflateThick() that deflate thin and thick triangles, respectively. The functions are mutually recursive since splitting a thin triangle produces a thin and a thick triangle, and splitting a thick triangle produces two thick and a thin triangle.
Listing C.24ch4 / Penrose
static final double phiInv = (Math.sqrt(5) - 1) / 2;
public static void deflateThin(
int depth, Point p, Point q, Point r,
List<List<Point>> triangles) {
if (depth == 0) {
triangles.add(List.of(r, p, q));
} else {
Point a = Point.lerp(p, q, phiInv);
deflateThin(depth - 1, r, a, q, triangles);
deflateThick(depth - 1, r, a, p, triangles);
}
}
public static void deflateThick(
int depth, Point p, Point q, Point r,
List<List<Point>> triangles) {
if (depth == 0) {
triangles.add(List.of(p, q, r));
} else {
Point a = Point.lerp(p, q, phiInv);
Point b = Point.lerp(p, r, phiInv);
deflateThick(depth - 1, b, a, p, triangles);
deflateThin(depth - 1, b, a, q, triangles);
deflateThick(depth - 1, r, b, q, triangles);
}
}
The tiling in Fig. 4.11 was constructed by deflating the two triangles of a thick rhomb seven times:
Point p = new Point(0, 0);
Point q = Point.polar(1, 72.0 * Math.PI / 180);
Point r = new Point(1, 0);
Point s = new Point(q.x() + r.x(), q.y() + r.y());
int depth = 1;
List<List<Point>> result = new ArrayList<>();
Penrose.deflateThick(depth, p, r, s, result);
Penrose.deflateThick(depth, p, q, s, result);
If desired, the resulting triangles can be recombined into rhombs.
Solution to Exercise 4.10
a)A natural way to model a parametric curve is an interface Curve with a method pointAt() that maps a number t to a Point. We define two additional methods tMin() and tMax(), which return the lower and upper limits of the parameter range. We can then implement the mystery curve as shown in Listing C.26.
Listing C.25ch4 / plotter
public interface Curve {
Point pointAt(double t);
double tMin();
double tMax();
}
Listing C.26ch4 / plotter
public class MysteryCurve implements Curve {
public Point pointAt(double t) {
return new Point(
cos(t) + cos(6 * t) / 2 - sin(-14 * t) / 3,
sin(t) + sin(6 * t) / 2 + cos(-14 * t) / 3);
}
public double tMin() {return 0;}
public double tMax() {return 2 * Math.PI;}
}
b)Figure C.6 shows the result of drawing both curves using 500 line segments each. For the mystery curve this yields a reasonably smooth curve, but since the butterfly curve is significantly more complex, at least twice as many line
segments would be necessary.
Figure C.6 Plots of the mystery curve and the butterfly curve using 500 line segments.
c)As an estimate for the straightness of this segment, we can use the distance between the midpointalong the curve \(\vec {C}((a+b)/2)\) and the midpoint of the two endpoints \((\vec {C}(a)+\vec {C}(b))/2\):
If the distance is less than a given threshold, we consider the segment straight; otherwise, we subdivide it further.
d)One possible program for subdividing arbitrary parametric curves is shown in Listing C.27. To draw a curve, the draw() method first
divides the parameter range between tMin() and tMax() into startIntervals pieces of equal size and then calls subdivide() to recursively subdivide each
interval. The subdivision process stops if the recursion depth has reached a given upper limit (maxSubdivisionDepth) or if the current segment is sufficiently straight, as measured by a parameter tolerance. The purpose of the maximum subdivision depth is to
prevent infinite recursion when plotting pathological curves that consistently fail the straightness test.
Listing C.27ch4 / plotter / CurvePlotter
// A class for drawing parametric curves using adaptive subdivision.
public class CurvePlotter {
public int startIntervals = 20;
public double tolerance = 1e-3;
public int maxSubdivisionDepth = 10;
public List<Point> draw(Curve curve) {
List<Point> points = new ArrayList<>();
double intervalSize =
(curve.tMax() - curve.tMin()) / startIntervals;
for (int i = 0; i < startIntervals; i++) {
double t0 = curve.tMin() + i * intervalSize;
double t1 = t0 + intervalSize;
subdivide(curve, t0, t1,
curve.pointAt(t0), curve.pointAt(t1),
0, points);
}
points.add(curve.pointAt(curve.tMax()));
return points;
}
private void subdivide(Curve curve, double tMin, double tMax,
Point pMin, Point pMax, int depth, List<Point> points) {
double tMid = (tMin + tMax) / 2;
Point mid = curve.pointAt(tMid);
double d = Point.distance(mid, Point.lerp(pMin, pMax, 0.5));
if (depth > maxSubdivisionDepth || d < tolerance) {
points.add(pMin);
} else {
subdivide(curve, tMin, tMid, pMin, mid, depth + 1, points);
subdivide(curve, tMid, tMax, mid, pMax, depth + 1, points);
}
}
}
e)The plots in Fig. 4.13 were created using the algorithm discussed in part (d) with a tolerance of 0.003. The same parameters were used for both curves.