Why is a raven like a writing desk?

In doing my modulation work with curves and ellipses, I extended the vecto function for drawing an ellipse to enable an oriented ellipse. Lately it occurred to me that this didn’t go far enough in terms of functionality, and I began wondering about how to draw part of an elliptical arc.

Vecto’s ellipse drawing function works by drawing the ellipse in quarters, with four cubic Bézier curves. Each curve’s control points are calculated according to the radii and kappa, the constant that is 4(√2 – 1)/3, as explained by G. Adam Stanislav. This works fine if you are drawing quarters, halves, or complete ellipses (or circles of course), but what I wanted was a function like this:

(defun approximate-elliptical-arc (cx cy a b theta eta1 eta2)
…)

With (cx, cy) being the centre, a and b being the radii, theta being the orientation, and eta1 and eta2 being the start and end angles specifying the arc. The return value would be a Bézier curve, or a list of curves, depending on the required accuracy.

So I did some research and found an interesting paper by Luc Maisonobe on approximating such elliptical arcs with lines, quadratic and cubic Bézier curves. The maths is straightforward, if rather heavy on the constants, and it was less than an hour’s work to code up the function. I now have functionality to draw arbitrary elliptical (or circular) arc segments.

The interesting thing is that this method of approximating ellipses is not the same as G. Adam Stanislav’s method: to put it another way, when I compute a quarter-circle arc with a single Bézier curve, the control points are not quite coincident with the kappa-method control points. The error in the circle approximation is still undetectable to the eye, but careful comparison of two large circles shows that the kappa-method one is slightly larger: the curves are slightly more “pushed out”.

G. Adam writes: “It is my contention that we get the closest to a real circle if we draw the one curve with properties already discussed, and with the center point of the curve lying at the same point the center point of the circumference of a true quadrant of a circle would lie.” I would be interested to know if he has a proof behind this contention, and/or an expression for the error in such a curve. The derivation of kappa is based around this idea that the midpoint of the curve should be coincident with the midpoint of the arc.

G. Adam’s method seems accurate. But the Maisonobe method can be used to subdivide the arc, obtaining successively closer approximations to the true arc by using multiple Bézier curves, and includes a method for calculating an upper bound on the error. This, together with the fact that it easily handles start and end angle specification of the arc, leads me to favour it.

Given Pythagorean triples that satisfy the Diophantine equation:

a² + b² = c²

where a, b and c share no common factors, one of a and b must be odd, the other must be even, and c is always odd.
First, note that squares of even numbers are always divisible by 4. [Lemma 1]

(2n)² = 4n²

And squares of odd numbers are of the form (4n + 1). [Lemma 2]

(2n + 1)² = 4n² + 4n + 1
= 4n(n+1) + 1

If a and b are both even, it is trivial to see that a, b and c share the factor 2.

c² = (2n)² + (2k)²
= 4n² + 4k²
= 4(n² + k²)

By Lemma 1, c is even. Therefore the solution reduces to the more primitive form:

n² + k² = (c/2)²

Consider a and b both odd:

c² = (2n + 1)² + (2k + 1)²
= 4n² + 4n + 1 + 4k² + 4k + 1
= 4n² + 4n + 4k² + 4k + 2
= 2(2(n² + n + k² + k) + 1)

This (by Lemmas 1 & 2) is not a perfect square, so this cannot be a solution. In fact, in this case c must be irrational:

Let p = n² + n + k² + k
c = √2√(2p + 1)

2p + 1 cannot have 2 as a factor. And since √2 is irrational, so is c. Therefore one of a and b must be odd and the other even.

c² = (2n + 1)² + (2k)²
= 4n² + 4n + 1 + 4k²
= 4n² + 4n + 4k² + 1
= 4(n² + n + k²) + 1

And c is odd (Lemma 2). QED.

Honestly, they’re incredibly interesting. Anyway, I’ll skip straight to the pièce de résistance:

This is a 4-curve cubic Bézier spline modulated onto an ellipse. The ellipse [a=4, b=3] is at an angle of π/4. C1 continuity of the complete curve is preserved. The flickr set tells the story of how I got here.

This is around 450 lines of my naive lisp, including class definitions and test code. So Gary, this is your g-code challenge!

The lisp code is object oriented (oh, and so much nicer than C++’s so-called object orientation). I rewrote the earlier code now that I knew what I was doing, and I added lines and polylines to the mix too (see the flickr set) so I can easily modulate whatever edge I want. You’ll notice if you look closely at the earlier attempts that they had a bit of a problem with c1 continuity, which is now fixed with the new code.

In closing, thanks, Zach!

After some code cleanup and generalisation, I can now modulate whole splines onto ellipses and onto splines themselves. Here is my simple 4-bezier spline modulated onto an ellipse:

And onto itself:

Repeatedly modulating a spline onto itself while varying the frequency parameter leads to some interesting and fractal patterns. Nice.

My efforts to equally subdivide a curve along its length have, in part, been leading to this. First, I extended the sampling to work with splines (made up of cubic Bézier curves with c1 continuity). This shot shows 4 curves put together to form a spline:

Next, I wrote some code to modulate a curve onto another curve. Here’s a curve, and the same curve with another curve modulated onto it.

Since I already had the binary search and interpolation code, it was just a matter of writing different samplers for ellipses and Bézier curves.

;; make a sampler function for a bezier curve
(defun make-bezier-sampler (p0 p1 p2 p3)
  (lambda (k)
    (decasteljau p0 p1 p2 p3 k)))

;; make a sampler function for an ellipse
(defun make-ellipse-sampler (center xradius yradius)
  (lambda (k)
    (let ((x (* xradius (cos k)))
          (y (* yradius (sin k))))
      (make-point (+ x (xcoord center)) (+ y (ycoord center))))))

Equal angle increments:

Equal circumferential distances:

I think I will use the same method as I do for Bézier curves to step along the circumference.

Another generalisation I had to make from the circle code is with respect to the normal. For a circle of radius r, centred on the origin, and parameterised by angle θ, a point on the circle is (r cos θ, r sin θ). And the normal is exactly the same: the point vector is normal to the circle.

For an ellipse, this is not the case. For an ellipse with half-dimensions a and b, centred on the origin, and parameterised by angle θ, a point on the ellipse is (a cos θ, b sin θ). But the normal is (b cos θ, a sin θ).

(Of course, the ellipse form degenerates to the circle form when a = b.)

Having conquered equidistant spacing along a Bézier curve, my thoughts now turn to the same problem for an ellipse. I have solved the problem for a circle of course, which is a special case of an ellipse. One would think that going from a circle to an ellipse would be mathematically easy: it’s easy to compute a point on an ellipse given centre and radii, and an ellipse is just a 2-way stretched circle, right?

Well, as with the Bézier curve, it’s not as simple as simply incrementing the angle parameter around the ellipse as one can with a circle, because obviously in that case the distance between successive points will vary.

And then one comes up against the rather startling fact that there is no simple exact equation for the circumference of an ellipse. No variant of 2πr here. As one reference puts it, “there are simple formulas but they are not exact, and there are exact formulas but they are not simple”. The exact formulas are infinite sums. The simple formulas can be glaringly inexact, depending on the ellipse.

I’m thinking on it some more. If, as I suspect, there proves to be no easy closed form solution to how much to vary the angle (or other alternative parameters of the curve) to achieve uniform spacing of radial points, I can fall back on the same solution as for the Bézier curve, i.e. sampling and interpolation.

As you know, a Bézier curve (here meaning specifically a cubic Bézier) is often used for drawing all kinds of things in vector graphics. It has the nice property that the endpoints and control points form a bounding box, and deCasteljau’s algorithm is a nice numerically stable way of evaluating the curve at a particular point. Of course, the same algorithm trivially gives the tangent to the curve at that point and from there it’s a hop and a skip to find the normal.

All well and good so far. The problem I wanted to tackle was that of spacing points equally along the curve (and incidentally, pushing these same points out along the normal), so that I can easily modulate a pattern (e.g. sawtooth, square wave) on to an arbitrary Bézier curve.

But there’s a slight problem with the naive method for doing this. Recall that a Bézier curve is in fact a parametric curve, specified by the parameter t which varies from 0 to 1 along the curve. The problem is that a linear increment of the parameter t does not correspond to a linear step along the curve length. For some curves, it’s not a bad approximation, but the more “curvy” the Bézier, the worse the discrepancy. This picture is a clear illustration of this problem.

The picture shows points (and points pushed along the normal) plotted on the curve, as t increments linearly from 0 to 1. As you can see, stepping along with constant increments of the parameter t gives the wrong result. These points are not equally spaced along the curve.

So far, I have a naive “correct” solution to this problem as follows: I sample the curve at a high frequency (say 1000). Each segment between points I then treat as a straight line, and I compute the length of the curve (simply the sum of the Pythagorean distance between successive points). I can then step along the curve linearly with respect to the curve length, and recover the parameter t (to actually evaluate the curve at a given point) by a simple calculation (binary search and interpolation between nearest neighbours). The result of this approach:

Much better.

Pascal’s Triangle (henceforth known as PT). You know – that thing you learned in maths. The Wikipedia entry, like most mathematical Wikipedia entries, reads (if your mathematical background is anything like mine) like “here’s a few things you might vaguely remember from school, oh and of course gleep = glorp”. Although I have to say, the PT entry is less opaque than most.

If you were lucky enough to have a “recreational maths” class then perhaps you studied PT more. I did have a recreational maths class in the 4th form (erm… 7th grade? I went to a public school with a nonstandard year numbering), but at the time I didn’t appreciate just how much higher level mathematics is intertwined (see trigonometry, complex numbers & calculus). Also, is it just me, or are curricula in general a bit light on actual number theory? I mean, we learn counting, and times tables, but after that it gets a bit fuzzy. Prime numbers? Some cursory investigation. But ISTM that students don’t really get to know the normal counting numbers very well. At least, I didn’t. I remember the occasions where we did take a lesson “off the curriculum” to study them as being some of my favourite maths lessons, e.g. the Hotel Cantor with a room for every positive integer, and the Infinity Bus Line which had a bus for every positive integer, each of which had a seat for every positive integer.

Anyway, PT. The thing I knew about but didn’t appreciate the significance of. Over the last week I’ve been tackling the Project Euler problems and there are a few concerning PT. (There are many inviting an exploration of number theory, which is also very cool). After solving 34%, I decided to skip down to problem 148, “Exploring Pascal’s Triangle”. Not many people have solved this, I thought: daunting! But also, this is my old friend PT. I’d like to get to know it better.

So after a couple of days of head scratching, I realised a way to attack the problem which would terminate before the heat death of the universe (always a useful approach to take). I also coded up some naive programs just to start exploring PT (as the problem says) – in the hope that they would give me some insight. Well, that and 3 or 4 pages of diagramming and scribbling later, I found out some quite interesting things. And I realised that I could solve the problem not just for multiples of 7, but for multiples of any prime. Now that’s cool.

I coded up a solution and ran it. It produced the wrong answer. I realised I’d missed something and corrected it. It produced the wrong answer. I tried again. Wrong answer. So I went to bed. Next morning, I realised my error, and also that I could simplify the program and speed it up. I coded it up again, and after ironing out a few non-algorithmic bugs, I ran it. Hooray! The right answer! I join the club of people who have solved problem 148.

Maths is cool.