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}
Improve this question

1 Answer 1

Reset to default
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

Improve this answer
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. Commented Dec 23, 2019 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. Commented Dec 23, 2019 at 19:51
  • Ah my bad, I miss understood the comment Commented Dec 23, 2019 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. Commented Dec 23, 2019 at 19:53
  • @makogan My comment wasn't clear, not your fault :) Commented Dec 23, 2019 at 19:55

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.