Archive for the ‘Lisp’ Category

And it was done

Saturday, March 1st, 2008

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:

ellipse-by-angle

Equal circumferential distances:

ellipse-by-circumference

More on ellipses

Saturday, March 1st, 2008

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.)

On ellipses

Saturday, March 1st, 2008

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.

Profiling with SBCL

Friday, February 29th, 2008

My Bézier curve code was not optimal, so I decided to learn how to profile with SBCL. In particular, I was not yet doing a binary search of the sampled points, neither was I doing an interpolation to recover the parameter t. Finally, I was sampling the curve at a fixed high frequency (probably too high for most cases). So I changed a list to an array, implemented binary search and interpolation to recover t, and set the sample frequency at 2x the number of division points.

SBCL profiling is very easy. The overall performance here is probably bound by the file write, which explains why the overall time for the top level function is more or less constant. But check out the other function timings and the “consed” column. The sample frequency change is clear in the number of calls to DECASTELJAU.

Before:

VECTO-TEST> (watch '(tenths-test-bezier "test.png"))
  seconds  |   consed   |  calls  |  sec/call  |  name  
----------------------------------------------------------
     0.278 | 26,952,592 |       1 |   0.278034 | TENTHS-TEST-BEZIER
     0.093 | 10,288,576 |  26,576 |   0.000004 | VEC-LENGTH
     0.037 | 10,236,064 |  26,525 |   0.000001 | POINT-DISTANCE
     0.015 |  1,607,920 |  12,624 |   0.000001 | INTERP-SINGLE
     0.003 |    385,024 |  32,990 |  0.0000001 | MAKE-POINT
     0.002 |          0 |  92,607 |  0.0000000 | XCOORD
     0.000 |          0 |       1 |   0.000000 | BEZIER-POINTS
     0.000 |  1,834,496 |      51 |   0.000000 | FIND-BEZIER-PARAM
     0.000 |     36,864 |       1 |   0.000000 | POLYLINE-LENGTH
     0.000 |     73,664 |       1 |   0.000000 | SAMPLE-BEZIER
     0.000 |      8,192 |      51 |   0.000000 | BEZIER-POINT-AND-NORMAL-AT
     0.000 |    339,968 |   1,001 |   0.000000 | DECASTELJAU
     0.000 |      8,192 |      51 |   0.000000 | NORMAL
     0.000 |     16,384 |      51 |   0.000000 | NORMALIZE
     0.000 |     16,384 |      51 |   0.000000 | SCALE
     0.000 |     40,960 |      51 |   0.000000 | ROTATE
     0.000 |  2,076,672 |   6,312 |   0.000000 | INTERP
     0.000 |     16,384 |      51 |   0.000000 | -POINT
     0.000 |          0 |  92,607 |   0.000000 | YCOORD
----------------------------------------------------------
     0.428 | 53,938,336 | 291,603 |            | Total

After:

VECTO-TEST> (watch '(tenths-test-bezier "test.png"))
  seconds  |   consed   |  calls |  sec/call  |  name  
---------------------------------------------------------
     0.439 | 26,773,984 |      1 |   0.438705 | TENTHS-TEST-BEZIER
     0.004 |    253,952 |    912 |   0.000005 | INTERP
     0.004 |     12,288 |    101 |   0.000035 | POINT-DISTANCE
     0.000 |     16,384 |      1 |   0.000000 | BEZIER-POINTS
     0.000 |     32,768 |    393 |   0.000000 | FIND-BEZIER-PARAM
     0.000 |     24,576 |      1 |   0.000000 | SAMPLE-BEZIER
     0.000 |     16,384 |     51 |   0.000000 | BEZIER-POINT-AND-NORMAL-AT
     0.000 |     32,768 |    101 |   0.000000 | DECASTELJAU
     0.000 |     57,296 |    152 |   0.000000 | VEC-LENGTH
     0.000 |          0 |     51 |   0.000000 | NORMAL
     0.000 |      8,192 |     51 |   0.000000 | NORMALIZE
     0.000 |     16,384 |     51 |   0.000000 | SCALE
     0.000 |          0 |     51 |   0.000000 | ROTATE
     0.000 |    175,856 |  1,824 |   0.000000 | INTERP-SINGLE
     0.000 |          0 |     51 |   0.000000 | -POINT
     0.000 |          0 |  2,535 |   0.000000 | YCOORD
     0.000 |          0 |  2,535 |   0.000000 | XCOORD
     0.000 |     16,384 |  1,166 |   0.000000 | MAKE-POINT
---------------------------------------------------------
     0.446 | 27,437,216 | 10,028 |            | Total

A couple of “it would be nice ifs” struck me while writing this code: IWBNI LOOP could collect into any sequence, not just a list; and IWBNI there was a CL equivalent of Haskell’s scanl. Both of these lacks could be solved (if they aren’t already and I’ve just overlooked them) but my CL macro skills aren’t all that yet.

Experimenting with (cubic) Bézier subdivision

Wednesday, February 27th, 2008

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.

bezier-subdiv-wrong

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:

bezier-subdiv

Much better.

Sneak preview

Tuesday, February 19th, 2008

This is a sneak peek at my current personal project. A bit of computerised heraldry.

8 way test

(“Sable a [ordinary] ermine.” Left to right, top to bottom: chief, fess, pale, cross, bend, bend sinister, chevron, saltire. Yeah, I know; charges on a bend should be bendwise. It’s a work in progress.)

Arc’s debut

Sunday, February 3rd, 2008

(Seems like this is just about what’s happened)

Paul Graham: Behold, Arc! (disclaimer: blah blah blah)

Rest of the Lisp-speaking world: Dude, weren’t you working for like 5 years on this? Now you give us a wrapper on Scheme?

PG: Look, I know it’s a bit pre-alpha, but give it a go! You might like it.

ROTLSW: Looks like a lisp with everything that makes actual real-world programming possible removed. Oh and some weird fairly pointless syntax changes?

PG: It made the code shorter. Look, a web app in five lines!

ROTLSW: …and five thousand lines of libraries. Even non-lispers can do that.

PG: *sigh* You’re just not getting it. All other lisps are Blubs now! I’m going back to John McCarthy’s papers…

ROTLSW: OK, call us when Arc can solve some real problems – see you in another 5 years then? Good luck with the web angle. We’ll be off playing with Vecto and other neat stuff in CL.

Lisping towards enlightenment

Monday, December 17th, 2007

Me: a C/C++ programmer of 12 years (professionally), a Lisp dabbler (programmer is too strong a word to use here) of, let’s say ~3 months. I’ve done enough functional stuff to be comfortable with that side, but I haven’t (hadn’t) yet written a project big enough to get me used to Lisp syntax.

The project: an LL(1) parser. That is, given a suitable grammar and set of legal lexemes, produce a function that will (top-down predictive) parse legal strings of said grammar. This isn’t rocket science by compiler standards, and I’m not sure it’s even likely that the grammar I have in mind to parse can even be made LL(1), but I have implemented precisely this in C++ before, so it’s a good test. Straight out of the Dragon Book section 4.4.

So far I have completed Lisp functions to produce a parsing table, using appropriate algorithms to construct First and Follow sets. This is about 110 lines of Lisp. Looking at the old C++ code (about 5 years ago) it comes to about 230 lines. From this I conclude that even a newb like me can gain immediate benefit from Lisp, and that I don’t yet know enough Lisp to make this about half the size it is (e.g. I’ve used no macros).

Time taken to get this far in Lisp: about 2 days. Time taken to get this far in C++: unknown, but the entire thing (grammar design, lexer, parser, bytecode compilation & VM execution) took on the order of a few months. The C++ version is production code that shipped in a real life game. Make of that what you will.

Even after 2 days, I’m already fine with the brackets. They really do not bother. The code is readable through indentation. I do sometimes run into problems with misplaced brackets when editing code, but it’s akin to a syntax error which is immediately found and fixed. I expect this to go away as I get more proficient.