NUBS and Polynomials, with graphic

As has happened before, while putting the finishing touches on the code for this post, I suddenly realized I needed to rewrite the code before I could be happy with it. So this post will be devoid of Perl 6 code, and just discuss what is going on in the (now-improved) graphic.

The basic idea I'm illustrating here is that a NUBS curve is basically a convenient way of merging a series of partial polynomial curves. In this graphic, the NUBS curve is in black. The polynomials are in red, green, and blue.



If you start where the red doesn't overlap the black, and then follow it to the black and follow the black beyond that, you can trace the curve and get an idea of what is going on. Basically, the black NUBS curve follows one section of the red polynomial, then smoothly switches to a section of the green polynomial, and finally switches to a section of the blue polynomial. It is piecewise polynomial, in other words, with smooth transitions. (It is possible to generate unsmooth transitions too with a NUBS, or even out-and-out discontinuities, but generally this is not wanted.)

The other great thing about NUBS is how natural they are to specify, because their control points conform roughly to the shape of the curve, and modifying them changes the curve in a fairly natural fashion. I don't have time to go into great depth here, but here's a quick comparison. First, here's the code to specify the curve in the picture:

my @control_points = (Vector.new(-1, -2),
                      Vector.new(1, 0),
                      Vector.new(1, 1),
                      Vector.new(0, 1),
                      Vector.new(1, 2),
                      Vector.new(1, 2),
                      Vector.new(1, 2));
my @knots = (-1, -1, -1, -1, 1, 2, 2, 3, 3, 3, 3);
my Nubs $nubs = Nubs.new(3, KnotVector.new(@knots), @control_points);

Then here's the red polynomial:

(0.194444444444444, 0.111111111111111) x^3 
+ (-0.916666666666667, -0.666666666666667) x^2 
+ (0.583333333333333, 1.33333333333333) x^1 
+ (0.694444444444444, 0.111111111111111) x^0

Actually, that's a fairly clean polynomial in this case, but even so, it's hard to work with. The bit we're interested in is parameterized from -1 to 1, so for instance, to find the starting point we need to evaluate the polynomial at x = -1:

- (0.194444444444444, 0.111111111111111)
+ (-0.916666666666667, -0.666666666666667) 
- (0.583333333333333, 1.33333333333333) 
+ (0.694444444444444, 0.111111111111111)
= (-1, -2)

Whereas the starting point of the NUBS version is just the first control point.

Ack. I could go on, but I fear I've exhausted the patience of the Perl readers out there. And I feel like I'm constantly butting up against the limits of Rakudo with this project, so I think I'll put it aside for a few months and tackle something else. Maybe try to whip up a grammar for ABC files, or something like that. Then return to Vector when after Rakudo has had a chance to mature a bit more...

Wiktory!

Rakudo has been fixed, and the code I've been trying to get work for a month works beautifully now! If I understand the fix properly, the problem was in Rat addition's call to Rat.new. The code was very dumb, so it always tried set the denominator of the new Rat to the product of the denominators of the two numbers being added. Rakudo's Ints are really Int32 right now (more or less), so if that product was equal to or greater than 2**31, it was autoconverted to Num. But the Rat.new which takes positional arguments takes Ints, so it wouldn't dispatch to that. Instead it would try to dispatch to the default autogenerated named argument form of Rat.new. But that only takes the implicit self parameter, and we were sending it self, Int, Num -- thus the "too many positional arguments: 3 passed, 1 expected" error!

Rakudo now autoconverts this case (424/61731 + 832/61731) to Num to avoid the overflow in the denominator. Obviously this is less than ideal -- a denominator of 61731 would work fine for this sum -- but it does work. And how! Gone are the mysterious crashes and errors. With the number of samples cranked up, the curves look beautiful.

Now I just need to do a bit of polishing to the output and figure out how to post it to the blog. I'm definitely feeling I've accomplished something cool in Perl 6....

PS Aha! Preview graphic, I'll explain what it means next time.

Aha!

A post-Thanksgiving dinner $work debugging session left me with some spare time to poke at my SVG.pm bug: that is, the "too many positional arguments: 3 passed, 1 expected" error I got when I tried to crank up the number of samples I was taking of the curve. I finally added enough says to track the error down to the simplest of operations: @N[$i] += $temp;, where $temp was 0.00686851014887172. This wasn't one of my fancy overloaded operators, it was a basic Perl 6 numeric operator. What could go wrong?

Well, I added a bunch more stuff to the say. And it turns out we are adding two Rats, 424/61731 and 832/61731. What could go wrong?

> say 424/61731 +  832/61731
too many positional arguments: 3 passed, 1 expected

If you look into infix:<+>(Rat, Rat), it just does a naive fractional add, relying on Rat.new to reduce the fraction. The problem here is that means the denominator of our new fraction starts its life as 61731 * 61731.

> say (61731 * 61731) div 246924
No applicable candidates found to dispatch to for 'infix:div'

I'm not sure what the best approach to fixing this is. Obviously I can change my code to do Nums. But I think the first step is to file a Rakudo bug.

Vector and SVG: Almost There

I noticed from #perl6 that chromatic fixed a memory leak in Parrot recently, and that fix made its way to Rakudo today. I thought, hey, I've been having a bus error. Could that be caused by an out-of-memory condition? So I upgraded Rakudo, and haven't seen the bus error since. (The "too many positional arguments: 3 passed, 1 expected" remain if I crank up the number of samples taken of each curve, but hopefully now I'll be able to track that down without being knocked around by bus errors.)

So, the first thing I did to get SVG working with Nubs and Polynomial was to write a simple class which converts from "normal" XY space to SVG coordinates (which I'm referring to as NM coordinates in this code).

class Vector { ... }
subset Vector2 of Vector where { $^v.Dim == 2 };

class SVGPad
{
    has Vector2 $.xy_min;
    has Vector2 $.xy_max;
    has Vector2 $.mn_min;
    has Vector2 $.mn_max;
    
    multi method new(Vector2 $xy_min, Vector2 $xy_max, Vector2 $mn_min, Vector2 $mn_max) 
    {
        self.bless(*, xy_min => $xy_min, xy_max => $xy_max, mn_min => $mn_min, mn_max => $mn_max);
    }
    
    multi method xy2mn(Vector2 $xy)
    {
        my $t = ($xy - $.xy_min).coordinates >>/<<  ($.xy_max - $.xy_min).coordinates;
        return $.mn_min + Vector.new(($.mn_max - $.mn_min).coordinates >>*<< $t);
    }

So the code to use this (for now) is just

class SVGPad { ... }
class Nubs { ... }
class Polynomial { ... }
sub MakePath($curve, Range $range, SVGPad $pad)
{
    my @points = RangeOfSize($range.from, $range.to, 10).map({$pad.xy2mn($curve.evaluate($_))});
    my $start = @points.shift;
    my $path = "M {$start.coordinates[0]} {$start.coordinates[1]}";
    for @points -> $v 
    {
        $path ~= " L {$v.coordinates[0]} {$v.coordinates[1]}";
    }
    return $path;
}

my @control_points = (Vector.new(-1, -2),
                      Vector.new(1, 0),
                      Vector.new(1, 1),
                      Vector.new(0, 1),
                      Vector.new(1, 2),
                      Vector.new(1, 2),
                      Vector.new(1, 2));
my @knots = (-1, -1, -1, -1, 1, 2, 2, 3, 3, 3, 3);
my Nubs $nubs = Nubs.new(3, KnotVector.new(@knots), @control_points);
my Polynomial $poly1 = $nubs.evaluate(0, Polynomial.new(0.0, 1.0));
my Polynomial $poly2 = $nubs.evaluate(1.5, Polynomial.new(0.0, 1.0));
my Polynomial $poly3 = $nubs.evaluate(2.5, Polynomial.new(0.0, 1.0));

my $pad = SVGPad.new(Vector2.new(-2.5, -2.5), Vector2.new(2.5, 2.5), 
                     Vector2.new(0, 0), Vector2.new(400, 400)); 
my $svg = svg => [
    :width(400), :height(400),
    path => [
        :d(MakePath($nubs, -1..3, $pad)), :stroke("blue"), :stroke-width(2), :fill("none")
    ],
    path => [
        :d(MakePath($poly1, -2..2, $pad)), :stroke("green"), :stroke-width(1), :fill("none")
    ],
    path => [
        :d(MakePath($poly2, 0..3, $pad)), :stroke("red"), :stroke-width(1), :fill("none")
    ],
    path => [
        :d(MakePath($poly3, 1..4, $pad)), :stroke("white"), :stroke-width(1), :fill("none")
    ],
];

say SVG.serialize($svg);

MakePath is a simple function which takes a "curve" (that is, an object which has an evaluate function which goes from t to x, y), a Perl 6 Range to evaluate it over, and an SVGPad, and returns a SVG path object. Then we set up a fairly simple NUBS curve, use the Nubs to Polynomial version of the evaluate function to generate the corresponding Polynomial for each segment of the curve, and output all four curves to SVG. This is still a crude first approximation, but it does work; if I knew how to include SVG in this post I could show it. Next time, hopefully.

Arrrgh!

Very briefly got out SVG output for a NUBS curve yesterday. Tried to make it a bit nicer looking, and had it promptly dissolve in a hail of bus errors. (And occasional "too many positional arguments: 3 passed, 1 expected" errors that I think must be incorrect.) It's now reminding me of my attempt to write an Euclidean geometry proof generator in Lisp back in '92 -- seemingly minor, innocent changes in the code lead to crazy changes in where it crashes.

Anyone have hints for how to approach debugging this sort of thing in Rakudo? I know it's just part of the frustration expected for an early Perl 6 adopter, but I'd really like to get this thing working...

NURBS Knot Vector Direction Done RIght

Last post I revealed how I'd gotten myself tangled up when trying to add a "do the right thing" mode to Nubs.evaluate. I'm going to repeat that code here, as I've added one of the Nubs.evaluate functions to the mix.

These functions became overly complex because I tried to overload the KnotBasisDirection enum with an extra layer of meaning. Pulling that extra layer out and rearranging things a tad gives you vastly better code:

Notice that we've eliminated an enum value here, two cases that needed to be checked for, and simplified a third line of code. The "trick" (almost too obvious to be called a trick) is to make the default value of $direction be determined by a called to Direction, rather than calling Direction with a special value that indicates it needs to do work.

It's the best of both worlds. By default, if you don't specify $direction it will do something reasonable. And for those rare weird cases where it is important (discontinuities in the curve), you can specify the direction to use for the evaluation.

Tomorrow: On the road to SVG in Perl 6.

Humbled by (Poor) Testing

I thought I was all ready to finish up with the Nubs / Polynomial thing last night. I started a blog post, then I looked at the code I wanted to show off. Well, it didn't seem clear enough to me to make me happy. Wait, I said to myself, this would make more sense as a Nubs member function. And that rewrite was so nice, I told myself it would be even better if I rewrote Nubs.Evaluate in the same fashion, because then the real beauty of this approach would shine out.

And then the excrement hit the fan, because after that perfectly straightforward rewrite, tests started failing left and right, with really weird errors. It took me about a half hour to figure out that this wasn't some strange Rakudo issue, it was a basic failure of my algorithm. For several weeks now, the code had been broken. But my tests were too clumsy to detect the failures until I stretched what I was doing really far out there.

Basically, the testing issue was this. The basis vector for the very first point of the NUBS curve should look like (1 0 0 0 ...). That means the very first point is just the very first control point. However, if in all your tests that first control point is some variant of (0, 0, 0, ...) (ie the origin), then a result of (0, 0, 0, ...) for the first point only shows you that the basis vector looks like (x 0 0 0 ...), where x can be anything at all. Whoops.

So it's back to the drawing board. Clearly I need more tests, and then I need to puzzle out where my algorithms have gone wrong. At least I'm pretty confident that I've got the Perl 6 issues worked out; now it's just a matter of fixing the math.

Polynomial & Vector: A Fly in the Ointment

This evening I was building my test case to make sure the Polynomial I am generating from the KnotVector basis actually yields the same results as the standard Nubs. Given an array of polynomials and an array of vectors, I wanted to multiply them with each other and sum the bunch. I wrote my $poly = [+] (@polys >>*<< @control_points);. When I ran it, I got back:

Ambiguous dispatch to multi 'infix:*'. Ambiguous candidates had signatures:
:(Polynomial $a, Any $b)
:(Any $a, Vector $b)

That's a beautiful error message, and completely appropriate. Both functions could do the job, each giving you their own answer. That is, (Polynomial $a, Any $b) where Any is a Vector would generate a Polynomial with Vector coefficients; (Any $a, Vector $b) where Any is a Polynomial would generate a Vector with Polynomial values. Both make sense, but for our purposes the former is probably best.

How to make that happen? The Polynomial and Vector classes are both completely independent of each other. Other than this potential intersection, they don't need to know about each other. After toying around a bit, I believe "is default" on the Polynomial multiply is the way to go.

Unfortunately, that gets me the error Method 'Num' not found for invocant of class 'Vector'. Possibly the problem is nested hyper multiply operators? Even when I map that operation, there are still two nested hyper multiply operators (both Polynomial and Vector multiply have them). Yet once again beautiful Perl 6 runs up against the awkward (and ever-improving) reality of Rakudo.

The Best Laid Schemes...

"The best laid schemes o' Mice an' Men, Gang aft agley" -- Robert Burns

I've hit on one of those programming conundrums where a quick diversion becomes a huge project. Let me explain what I'm going for; it will probably take three or four blog posts to cover all the details.

My basic idea for a cool Perl 6 twist on knot vectors was evaluating the basis function to create a polynomial rather than a set value. The basis function defines a piecewise polynomial, with a continuous polynomial between each pair of distinct neighboring knots. So if the knot vector is (0, 0, 0, 1, 2, 3, 3, 3), that defines three polynomials: from 0 to 1, from 1 to 2, and from 2 to 3.

For any reasonably complicated knot vector, those polynomials are pretty tricky to calculate by hand. The NURBS Book approach is a change of basis, which works well but is not particularly mathematically enlightening. But the knot vector basis calculation function N we've defined can already do this: if it is passed a reasonably capable polynomial object for the $u parameter. This is the payoff for not forcing $u to be a Num; we can pass a polynomial instead and get a polynomial out.

Well, almost. We have two KnotVector.N functions, one for the $p == 0 case and one for $p > 0. The latter will work perfectly when passed a polynomial for $u. (At least in theory; there may be issues with Rakudo bugs, of course.) The former, however, actually needs to know what the numeric value of $u is. Originally I planned to somehow overload the polynomial class to have a value for this purpose. But as far as I know, to make that work, I'd have to overload prefix:<+> -- and last time I checked, that didn't work in Rakudo.

You could fix that by changing the second N to take an impulse array, instead of calculating it using the value of $u. But once you're thinking about doing that, it's time to face the facts: even the hypothetical "cool" implementation of N we have is grotesquely inefficient for reasonably complicated knot vectors.



To calculate the basis efficiently, you should only calculate the values that have a chance of being non-zero. And you should get the impulse value by doing something like a binary search.

So before I can get where I want to go with knot vectors, I need to have a working polynomial class and a working binary search. Plus I think it might be worthwhile to have a short post discussing the .Str and .perl functions for KnotVector and Nubs...

NURBS Knot Vectors In Perl 6 (part 3)

Well, instead of being a slick demonstration of how powerful Perl 6 is, it has turned into a slog to find something the current Rakudo is happy with. The beautiful hyper-equation of my last post has turned into a dull and confusing series of little hyperoperations (and one which got expanded because the hyperoperator just wouldn't work).


I think all the bugs I've encountered here have already been reported, but I'm not sure. If anyone has suggestions on how to ably search RT, I'm all ears. My attempts so far have been painfully clumsy, to the point where my inclination is either to not report if I think it's likely to have been reported so far, or just report it without searching.

At any rate, let me finish by quickly sketching out a NUBS (Non-Uniform BSpline) curve class. This isn't quite as powerful a tool as a NURBS (Non-Uniform Rational BSpline), but it's still pretty powerful, and it's dead easy to implement the basic evaluate function with the tools we have so far.


Combined with the Vector class, this allows you to easily define N-dimensional NUBS curves. But we don't force the control points to be Vectors -- you can get a NUBS representation of polynomials over reals or complexes easily using this too, and if you have your own Vector class you'd prefer, that should work too, as long as you have scalar multiplication and vector addition implemented using * and +.

I do have one interesting twist left for another post on the KnotVector class, something I've always wanted to do that wasn't very practical in C++. And then I will work on expanding Nubs, and implmenting a proper Nurbs class, and looking at more complicated structures for each (surfaces, etc).

NURBS Knot Vectors In Perl 6 (part 1)

I'm not going to provide a lot of background for the wheres and whys of Non-Uniform Rational B-Spline (NURBS) knot vectors yet. Basically, a NURBS object is one or more knot vectors which define the parametrization and polynomials of the object, along with a bunch of Vectors to define the frame of the shape. I'm planning on spending a couple of posts working on getting a good knot vector class built up.

Confusingly enough, knot vectors have nothing directly to do with the sort of Vector we have been working with so far. From, say, a Perl 6 perspective, it would probably be more correct to call them knot arrays or knot lists. But knot vector is the traditional term, so that is what we will use.

The knot vector is a series of numbers in non-descending order, some of them repeated, that define the parameter range of the NURBS object we are looking at, and where it does interesting things. That is to say, adding a knot to the middle of the knot vector essentially splits one polynomial span the NURBS object is defined over into two joined polynomials. It's a fairly simple basis for a lot of powerful geometry.

Instead of starting by defining a knot vector class which can be the basis for NURBS objects, we're going to start by implementing the definition of the B-spline basis functions directly. (Note that the Wikipedia version there specifies it using different variable names -- our notation comes from The NURBS Book.) This would be horribly inefficient for real-world calculations, but we're just going to use it to establish good tests for our real class when that is ready.

Here's a straight forward implementation of the basis function formulas.



There are two tricky bits here. The first is that the standard definition can result in dividing zero by zero, which is defined to result in zero in this case. To keep the function implementation clean and simple, we define a "O/" operator which implements this exception to the standard rules of division.

The second wrinkle is that the standard definition has one curious blind spot. The equation does not work out when you reach the very last point! That is to say, the B-spline basis defined the standard way on a knot vector whose first value is A and last value is B covers the range [A, B). (Actually, I'm using first and last rather loosely there -- that's only true with the most common kind of knot vectors.) We define a KnotBasisDirection argument to N allows you to switch to the symmetrical definition, which is defined over (A, B]. The two definitions are identical over the interior range (A, B). The KnotBasisDirection enum has two values, Left and Right, corresponding to these choices. That's the first time I've defined an enum in Perl 6, but it seems straightforward enough.