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.

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.

Polynomial & Vector: Together At Last

After a weekend spent studiously ignoring Perl 6, I've finally dug back in and sorted out the problems in the KnotVector / Nubs implementation. At least, every test is back to passing, and I've added a number of harder tests which pass as well. The good news is I like the code I have now a lot better than the code I thought I was ready to post last time around. (Still needs some work, though, IMO.)

Without further ado, here's the heart of what I've been trying to get at.

multi method Evaluate($t, KnotBasisDirection $direction = Left)
{
    my $n0 = $.knot_vector.N0_index($.degree, $t, $direction);
    return [+] ($.knot_vector.N_local($n0, $.degree, $t) 
                >>*<< @.control_points[$n0 .. ($n0 + $.degree)]);
}

multi method Evaluate($base_t, $actual_t, KnotBasisDirection $direction = Left)
{
    my $n0 = $.knot_vector.N0_index($.degree, $base_t, $direction);
    return [+] ($.knot_vector.N_local($n0, $.degree, $actual_t) 
                >>*<< @.control_points[$n0 .. ($n0 + $.degree)]);
}

Notice these two routines are exactly the same, except where the second uses $base_t and $actual_t, the first uses just $t both times. This allows us to always set $base_t with an actual number (so that it can be used to calculate $n0), but pass anything that logically works for $actual_t. If $actual_t is a number, then Evaluate will calculate the curve's value at that parameter value. If it is a Polynomial representing a single variable (call it what you will, I keep switching between u, t, and x in my thinking), then Evaluate will return the Polynomial representing the section of the curve specified by $base_t. If (as I keep toying around with), you passed a class representing a closure you can do math on, then you'll get a fancy closure that can evaluate that section of the curve. And so on and so forth. It's not tied to a particular implementation of Polynomial; all that's required is something that can handle a few basic math functions.

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: Debugging Update

I've tracked down the next bug mentioned at the end of my last post, and again it comes as a side effect of the hack to make it easier to code Polynomial addition. In that hack, we extend the Polynomial coefficients array with 0. But if the coefficients are Vectors, the actual extension value should be Vector.new(0, 0, 0). Just plain 0 means we eventually try to add 0 to a Vector, and since Vector doesn't support addition with a scalar, it falls back to a normal addition operator, which tries to convert Vector to a Num.

Thoughts: I don't see any obvious way to figure out the correct "zero" to extend the coefficients with. Nor do I see any obvious way to allow you to add 0 to a Vector. (At least without opening up general vector plus scalar math, which is a bad idea IMO.)

It also suggests that it may be a very good idea to go beyond implementing the operators that work to implement dying versions of the operators that shouldn't be allowed, because allowing Perl 6 to fall back to its own operators will produce less useful error messages. (Or worse, they might even "work".)

I can see one more potential error like this lurking in the Polynomial.prefix<-> code...

Polynomial & Vector: Debugging

So, I was just going to post on working around the bug mentioned in my last post, when I noticed Moritz++ had commented with a suggested debugging approach. And we're off and running!

multi method Num()
{
die "Cannot call Num on Vector!";
}

It's not the most elegant error message, but it will do. Adding that switches from the fairly useless Method 'Num' not found for invocant of class 'Vector'
in Main (file src/gen_setting.pm, line 206) message to the extremely interesting

Cannot call Num on Vector!
in method Vector::Num (file lib/Vector.pm, line 24)
called from method Polynomial::new (file , line )
called from sub infix:* (file lib/Polynomial.pm, line 106)
called from Main (file , line )

Aha! Not a Rakudo bug at all, this is a legit issue. When I added the fix to remove trailing zero coefficients, I wrote @x[*-1] == 0 to check for them. Odds seem very good that == calls .Num. So... a change to .abs here is probably more realistic, anyway. (Maybe? I need to ponder that.)

Hmmm... that moves the error elsewhere, to (according to the backtrace) infix:+ (file , line ). I need to go to bed now, so the final debugging will have to wait. But at a minimum, Moritz's suggestion of adding a .Num that dies appears to be a valuable Perl 6 debugging trick.

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.

Polynomial and <<+>>

Defining a one-variable polynomial class is actually a pretty straightforward thing. We use a standard Perl 6 array to hold the polynomial's coefficients, with the Nth element of the array representing the coefficient of x^N. You can check out the complete source if you'd like; for the most part, Perl 6's nifty hyperoperators and list operations make the implementation simple.

One example where that isn't true is simple polynomial addition. In principle, this is the simplest thing possible; you just do a pairwise add on all the coefficients of the two polynomials. The only catch is when the polynomials have different degrees, so the arrays have different sizes. In that case we'd like to pad the array with zeros. Unfortunately, <<+>> pads with the last element of the shorter array, with is completely mathematically unsuitable in this case. So here was my first stab at implementing Polynomial addition:


That sort of thing works, but seems more like it belongs in the ugly world of C++ than the elegant world of Perl 6. The subject came up on the #perl6 channel yesterday, and pmichaud++ came up with a nice workaround: $a.coefficients, 0 <<+>> $b.coefficients, 0. As it is, this adds an additional zero coefficient each time it is called, which in the long run will be trouble. But it can be counteracted by adjusting the Polynomial constructor to remove "trailing" zero coefficients. (This is how Polynomial.pm is currently implemented.) (Errr... perhaps I spoke too soon; I'm having difficulty successfully deleting the trailing zeros. Sigh.)

On the IRC chat, colomon suggested using an adverb to tell <<+>> how to extend. His suggestion of :extend-with seems like a weak name to me, but I do think it seems like the proper Perl 6ish way to handle this issue.

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