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

Vector: Quick update

I meant to post on my polynomial trick, which now works, but instead got into a long run of small changes this evening. So let me just say that abs has been moved to the setting, so tonight I wrote Vector.abs, which makes Test.pm's is_approx work to compare two Vectors. It feels like Perl 6 is becoming more useful every day.

I also put in a hack workaround for the rakudobug (reported this evening) turned up in t/05-search.t, which means all the Vector tests pass again.

I'm also starting to ponder using SVG to output some NURBS (well, NUBS for now) curves to provide a graphic illustration of what my new code here is doing...