Skip to main content

splines - How to make the result of Piecewise be a closed interval?


Description


This question comes from two questions. Namely Q1 and Q2


The defintion of B-Spline basis function as shown below:


Let $\vec{U}=\{u_0,u_1,\ldots,u_m\}$ a nondecreasing sequence of real numbers,i.e, $u_i\leq u_{i+1}\quad i=0,1,2\ldots m-1$


$$N_{i,0}(u)= \begin{cases} 1 & u_i\leq u

  NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),u_] /;
i <= Length[knots] - 2 :=

Piecewise[
{{1, knots[[i + 1]] <= u < knots[[i + 2]]}},0]

coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])

NBSpline[i_Integer, p_Integer, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),
u_] /;p > 0 && i + p <= Length[knots] - 2 :=
Module[{init, res},
init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];

res = First@
Nest[
Dot @@@
(Thread@
{Partition[#, 2, 1],
With[{m = p + 2 - Length@#},
Table[
{coeff[u, k + 1, k + m + 1, knots],
coeff[u, k + m + 2, k + 2, knots]}, {k, i, i + Length@# - 2}]]}) &,
init, p]

]

Compare to built-in function BSplineBasis


 sortResult[x_ /; x == 0] := 0;
sortResult[res_] := MapAt[SortBy[#, Last] &, res, 1]

pts = {{0, 0}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};

Comparsion 1



 NBSpline[#, 0, knots, u] & /@ Range[0, 8] // Simplify


enter image description here



 PiecewiseExpand@BSplineBasis[{0, knots}, #, u] & /@ Range[0, 8]


enter image description here




Comparsion 2


sortResult /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7] // Simplify)


enter image description here



sortResult /@ (PiecewiseExpand@BSplineBasis[{1, knots}, #, u] & /@ 
Range[0, 7] // Simplify)



enter image description here





By the comparsion, I found that the result of BSplineBasis is always a closed interval, whereas the result my function NBSpline is a half closed- half open interval.


My trail:


postProcess[res_ /; res == 0] := 0;

postProcess[res_] :=
Module[{interval, pos, expr},
interval =

First@Simplify@res;
pos =
Position[
SortBy[interval, Last], Less];
expr =
ReplacePart[
SortBy[interval, Last], pos[[-1]] -> LessEqual];
Piecewise[expr, 0]
]


Test


 postProcess /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7])


enter image description here



Question:


Is there any good/simple method(strategy) to deal with this problem? I think my method is complex.



enter image description here





Answer



a fix, specially treating the end interval:


 NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),
u_] /; i <= Length[knots] - 2 := (
If[ knots[[i + 2]] != knots[[-1]],
Piecewise[{{1, knots[[i + 1]] <= u < knots[[i + 2]]}}, 0],
Piecewise[{{1, knots[[i + 1]] <= u <= knots[[i + 2]]}}, 0] ]);
coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])

NBSpline[i_Integer, p_Integer,
knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &), u_] /;
p > 0 && i + p <= Length[knots] - 2 :=
Module[{init, res},
init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];
res = First@
Nest[Dot @@@ (Thread@{Partition[#, 2, 1],
With[{m = p + 2 - Length@#},
Table[{coeff[u, k + 1, k + m + 1, knots],
coeff[u, k + m + 2, k + 2, knots]}, {k, i,

i + Length@# - 2}]]}) &, init, p]]
pts = {{1/10, 2/10}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};
BSplinePlot1[pts : {{_, _} ..}, knots_, opts : OptionsPattern[Plot]] :=
Module[{p = Length@First@Split[knots] - 1, a, b},
{a, b} = {First[knots], Last[knots]};
f[u_] =
Evaluate@
Simplify@
Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &,

pts] /. x_Piecewise :> Piecewise[x[[1]], Null];
ParametricPlot[f[u], {u, a, b}, opts, PlotPoints -> 1000]]
BSplinePlot1[pts, knots, ImageSize -> 600]

this bit: /. x_Piecewise :> Piecewise[x[[1]], Null] fixes the catchall behavior of Piecewise so you will get an error if you try to evaluate out of bounds. (not necessary for this example )


here is f[u] , note that we have closed just the end interval..


enter image description here


note as a practical matter @belisarius approach works as well, so long as you apply at the last stage:


 f[u_] = Evaluate@
Simplify@

Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &,
pts] /. Less :> LessEqual /. x_Piecewise :> Piecewise[x[[1]], Null]

enter image description here


since the function is continuous it doesn't really matter if you close all the internal intervals.


Comments

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...

list manipulation - Selecting multiple columns from a matrix?

Sample data: data = { {{2013, 1, 1}, 24.13, 167.67, 231.82}, {{2013, 1, 2}, 32.15, 170.92, 225.99}, {{2013, 1, 3}, 35.43, 172.68, 221.67}, {{2013, 1, 4}, 36.73, 173.05, 218.32}, {{2013, 1, 5}, 58.19, 165.96, 197.05}, {{2013, 1, 6}, 69.99, 163.50, 187.52}, {{2013, 1, 7}, 71.37, 154.21, 175.58}, {{2013, 1, 8}, 72.51, 149.66, 163.25}}; I want a DateListPlot with three graphs, so for a matrix formed by columns 1 and 2, one for columns 1 and 3, and 1 for columns 1 and 4. At the moment I'm using this code: data2 = Transpose[{data[[All, 1]], data[[All, 2]]}]; data3 = Transpose[{data[[All, 1]], data[[All, 3]]}]; data4 = Transpose[{data[[All, 1]], data[[All, 4]]}]; DateListPlot[{data2, data3, data4}, Joined -> True, Filling -> {3 -> {1}}] but I have a hunch that this can be done more efficiently. I don't like the Transpose s in particular. Any ideas? edit (for extra credit) What if I need to multiply the second column by 2, which in my solution is simp...