2

I am currently trying to make a general B spline basis plotter with meta post.

The definition of a basis function is on wikipedia under properties.

enter image description here

The knots of a b spline should have the pattern 0's a number of times equal to the order of the spline, ascending values for a while, last value repeated to match the order

So for example for a spline of order 3: [0,0,0,1,2,3,4,5,6,6,6] makes a knot vector.

To this effect I tried coding that exact formula and printed the knot vector I fabricate to make sure I did things properly:

enter image description here

The above is one of the order 2 splines, it should look like a liner rising and then decreasing function (a pyramid) but instead it looks like an order 1 basis. I tried changing the order, but all of my functions are returning the same pattern.

Like, for example an order 3 function: enter image description here

My code is as follows;

\documentclass[border=10cm]{standalone}

\usepackage{luamplib}
\mplibnumbersystem{double}
\usepackage[margin=0.5cm]{geometry}

\begin{document}
{\centering
\begin{mplibcode}
u:=1cm;

numeric knots[];

vardef make_knots(expr order, num)= 

    p_num = num - 1;
    v = 0;
    for i=0 upto order-1:
        knots[i] := v;
    endfor;
    for i=0 upto p_num - order:
        v := v+1;
        knots[i + order] := v;
    endfor;
    for i=0 upto order-1:
        knots[i + p_num + 1] := v;
    endfor;

    for i=0 upto order + p_num - 1:
        tmp := round(knots[i]);
        label.top(textext("\huge$"& decimal(tmp) &"$"), (2*i*u,6*u));
    endfor;
enddef;

vardef calculate_basis(expr t, i, order)= 
    numeric ret;
    if order=1:
        if (t >= knots[i]) and (t < knots[i+1]):
            ret = 1;
        else:
            ret = 0;
        fi;
    else:
        ret = 
        ((t-knots[i]) / (knots[i] + order - knots[i])) * calculate_basis(t, i, order-1) 
        + 
        ((knots[i + order + 1] - t) / (knots[i + order + 1] - knots[i+1])) * calculate_basis(t, i+1, order-1);
    fi;
    ret
enddef;

color darkred, darkyellow, darkgreen, lightblue, brown, pink, orange; 
darkred := (0.8,0.0,0.0);
darkyellow := (1.0,0.8,0.0);
darkgreen := (0.0,0.6,0.0);
lightblue := (0.0,0.8,1.0);
brown := (0.5, 0.1, 0.1);
pink := (1, 0.0, 0.8);
orange := (1, 0.4, 0.0);
%ignore first parameter while denbugging
vardef plot_basis(expr j, order, color)= 
    res := 100;

    for i=0 upto res-1:
        fraction := (i/res) * 6; %multiply by biggest knot value
        valS := calculate_basis(fraction, 1, 3);

        save pointS;
        pair pointS;
        pointS := ((i) * 2 / res, valS) * u * 5; %scale plot to make it visible
        draw pointS withpen pencircle scaled 5bp withcolor color;

    endfor;
enddef;

% Start figure
beginfig(0);

color colors[];
colors[0] = darkred;
colors[1] = darkyellow;
colors[2] = darkgreen;
colors[3] = lightblue;
colors[4] = pink;
colors[5] = brown;
colors[6] = orange;

make_knots(2, 8);
plot_basis(2, 2, colors[0]);

endfig;
\end{mplibcode}
\par}
\end{document}
2

I can't quite explain why, but I do have a fix (I think) and maybe also a couple of fixed typos. In your recursion (see comments):

vardef calculate_basis(expr t, i, order)=
    numeric ret;
    if order=0: % I think this should be zero to match your screenshot
        if (t >= knots[i]) and (t < knots[i+1]):
            ret = 1;
        else:
            ret = 0;
        fi;
    else:
        ret = 
        ((t-knots[i]) / (knots[i + order] - knots[i])) * calculate_basis(t, i, order-1) % fixed typo in denominator of first fraction
        + 
        ((knots[i + order + 1] - t) / (knots[i + order + 1] - knots[i+1])) * calculate_basis(t, i+1, order-1);
    fi;
    ret
enddef;

In the else case, I think that ret is somehow taking the value of the final order 0 case evaluated in the recursion. I don't know anything about b-splines, and not much more about how the grouping, recursion calls, and equation solving interact here. Why the recursion works is confusing to me because it looks like it expands to something like

ret=stuff1*(ret=0)+stuff2*(ret=1)

which doesn't make any sense to me. Hopefully someone that can explain the why comes along and does so.

Either way, this is how I would have written the recursion:

vardef calculate_basis(expr t, i, order)=
    if order=0:
        if (t >= knots[i]) and (t < knots[i+1]):
                    1
        else:
                    0
        fi
    else:
        ((t-knots[i]) / (knots[i + order] - knots[i])) * calculate_basis(t, i, order-1) % typo
        + 
        ((knots[i + order + 1] - t) / (knots[i + order + 1] - knots[i+1])) * calculate_basis(t, i+1, order-1)
    fi
enddef;

and after making a final change

    valS := calculate_basis(fraction, 3, 1); % changed to first order

the following seems to produce something closer to what you want (I'm not sure where the remainder of the fixed values in the plot basis loop come from). Either way, hopefully this will give you a start:

\documentclass[border=10cm]{standalone}

\usepackage{luamplib}
\mplibnumbersystem{double}
\mplibtextextlabel{enable}
\usepackage[margin=0.5cm]{geometry}

\begin{document}
{\centering
\begin{mplibcode}
u:=1cm;

numeric knots[];

vardef make_knots(expr order, num)= 
    p_num = num - 1;
    v = 0;
    for i=0 upto order-1:
        knots[i] := v;
    endfor;
    for i=0 upto p_num - order:
        v := v+1;
        knots[i + order] := v;
    endfor;
    for i=0 upto order-1:
        knots[i + p_num + 1] := v;
    endfor;

    for i=0 upto order + p_num - 1:
        label.top("\huge$"& decimal(round(knots[i])) &"$", (2*i*u,6*u));
    endfor;
enddef;

vardef calculate_basis(expr t, i, order)=
    if order=0:
        if (t >= knots[i]) and (t < knots[i+1]):
                    1
        else:
                    0
        fi
    else:
        ((t-knots[i]) / (knots[i + order] - knots[i])) * calculate_basis(t, i, order-1) % typo
        + 
        ((knots[i + order + 1] - t) / (knots[i + order + 1] - knots[i+1])) * calculate_basis(t, i+1, order-1)
    fi
enddef;

color darkred, darkyellow, darkgreen, lightblue, brown, pink, orange; 
darkred := (0.8,0.0,0.0);
darkyellow := (1.0,0.8,0.0);
darkgreen := (0.0,0.6,0.0);
lightblue := (0.0,0.8,1.0);
brown := (0.5, 0.1, 0.1);
pink := (1, 0.0, 0.8);
orange := (1, 0.4, 0.0);

vardef plot_basis(expr order, color)= 
    res := 100;
        pair pointS;
    for i=0 upto res-1:
        fraction := (i/res) * 6; 
        valS := calculate_basis(fraction, 3, 1); % changed to first order
        pointS := ((i) * 2 / res, valS) * u * 5; 
        draw pointS withpen pencircle scaled 5bp withcolor color;
    endfor;
enddef;

% Start figure
beginfig(0);

color colors[];
colors[0] = darkred;
colors[1] = darkyellow;
colors[2] = darkgreen;
colors[3] = lightblue;
colors[4] = pink;
colors[5] = brown;
colors[6] = orange;
make_knots(2, 8);
plot_basis(1, colors[0]);
endfig;
\end{mplibcode}
\par}
\end{document}

enter image description here

7
  • Thanks a lot, I will try it when I get home. As to the explanation, each spline of order 0 is 1 on a given interval and 0 every where else, so the combinations are actually (stuff * 0) + (stuff * 0) or (stuff * 1) + (stuff * 0) or (stuff * 0) + (stuff * 1). Hopefully that helps with the confusion. If you see on your result (which is what I expect to see) the graph is first 0, then raises to 1 on an interval, then decreases to 0 again on the next interval. – Makogan Dec 23 '19 at 16:14
  • 1
    @Makogan I sort of understood what the output was supposed to be; the confusion was more on the metapost side of things. In metapost x=2 is an equation that metapost solves for x, whereas x:=2 is an assignment. In the loop stuff*0 isn't what's happening, it's stuff*(ret=0) and the ret=0 part doesn't return a value as far as I'm aware so I don't know how metapost is evaluating what's happening in the else part of the main conditional other than it's somehow assigning to ret one of the values from the base case. – Scott H. Dec 23 '19 at 19:51
  • Ah my bad, I miss understood the comment – Makogan Dec 23 '19 at 19:53
  • @makogan Another way to get the correct result is using assignment in the conditional in the final else: ret := ((t-knots[i]) / (knots[i] + order - knots[i])) *... or assignments everywhere. – Scott H. Dec 23 '19 at 19:53
  • @makogan My comment wasn't clear, not your fault :) – Scott H. Dec 23 '19 at 19:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.