1

I'm trying to implement the recursive definition for B-Splines in c# but I can't get it right. Here's what I've done:

public static Double RecursiveBSpline(int i, Double[] t, int order, Double x)
{
    Double result = 0;

    if (order == 0)
    {
        if (t[i] <= x && x < t[i + 1])
        {
            result = 1;
        }
        else
        {
            result = 0;
        }
    }
    else
    {
        Double denom1, denom2, num1, num2;

        denom1 = t[i + order + 1] - t[i + 1];
        denom2 = t[i + order] - t[i];

        if (denom1 == 0)
        {
            num1 = 0;
        }
        else
        {
            num1 = t[i + order + 1] - x / denom1;
        }

        if (denom2 == 0)
        {
            num2 = 0;
        }
        else
        {
            num2 = x - t[i] / denom2;
        }

        result = num1 * RecursiveBSpline(i + 1, t, order - 1, x) 
            + num2 * RecursiveBSpline(i, t, order - 1, x);
    }

    return result;
}

And here is how I call the function:

Double[] vect = new Double[] { 0, 1, 2, 3 };

MessageBox.Show(BSpline.RecursiveBSpline(0,vect,2,0.5).ToString());

I should see 0,125 on the screen, instead I get 0,25. The two denominator variables are used to check if they equal 0 and if they do, the number should be set to 0 by definition. Can someone point out where I'm getting this wrong?

Lee Taylor
  • 7,761
  • 16
  • 33
  • 49
Gábor Birkás
  • 623
  • 3
  • 13
  • 20
  • Floating-point equality comparisons (`denom1 == 0`, `denom2 == 0`) [can be a little tricky](http://floating-point-gui.de/), because of e.g. slight rounding errors. Even when you *think* a number ought to be 0, can you be sure that no arithmetic operation (e.g. the minus division operator) introduced a minute rounding error? – stakx - no longer contributing Oct 19 '13 at 15:27

1 Answers1

1

Bear in mind, that the mathematical and logical operators in C# have a precedence order. Your second solution works fine if you put the right terms in braces (explanation follows). This line:

num2 = x - t[i] / denom2;

should be changed to:

num2 = (x - t[i]) / denom2;

and so on. Then the result is as desired: 0.125

The division operator has a higher order precedence as the addition operator. To affect the execution order use braces (everything in braces will be evaluated at first):

var r1 = 2 + 2 / 2; // Step1: 2 / 2 = 1 Step2: 2 + 1 Output: 3
var r2 = (2 + 2) / 2;   // Step1: (2 + 2) = 4 Step2: 4 / 2 = 2 Output: 2
keenthinker
  • 7,645
  • 2
  • 35
  • 45
  • Thank you very much, I've never thought of that! I've edited my original question in the meantime but your solution works well! – Gábor Birkás Oct 19 '13 at 15:13
  • Maybe it's not sorted out after all :D I used your solution to the recursive function above and I get correct results with 0<=x<=1. But if x>1 I should get 0 but this gives me 0.75 for x=1.5 and i=0. – Gábor Birkás Oct 19 '13 at 15:38
  • See the comment from stakx. Maybe the conditions aren't always met correctly. Is x > 1 an allowed input? – keenthinker Oct 19 '13 at 15:45
  • I've read that comment and that shouldn't be the problem. And x>1 is of course an allowed input but the recursive definition should produce 0. Maybe I'm mixing up some indexes but for the life of me I can't figure out what's wrong and I can't find an implementation of this in c# despite that it's a really famous definition. – Gábor Birkás Oct 19 '13 at 15:52
  • Have you evaluated the result manually for x > 1? From where did you took the definition of the function? Some of your terms / indices are not the same as for example the definition given in the english wikipedia article for a B-Spline. – keenthinker Oct 19 '13 at 16:02
  • That's because wikipedia's definition worked even worse. Those indexes are incorrect in a zero based context. I took the definition from here: [NURBS](http://web.cs.wpi.edu/~matt/courses/cs563/talks/nurbs.html) – Gábor Birkás Oct 19 '13 at 16:05
  • Everything's fine. What I did was I misinterpreted the meaning of symbols in wikipedia. Thanks again for all the help! – Gábor Birkás Oct 19 '13 at 16:30