B-spline
- Task
Generate a B-spline curve with a list of 12 points and plot or save image.
Coordinates of control points:
start=171,171 1 185,111, 2 202,109, 3 202,189 4 328,160 5 208,254 6 241,330 7 164,252 8 69,278 9 139,208 10 72,148 end=168,172
Rules!!!!
Do not use third party libraries or functions
- See also
Julia
Choose BSpline D of 2, ie degree 1. <lang julia>using Graphics, Plots
Point(t::Tuple) = Vec2(Float64(t[1]), Float64(t[2])) const controlpoints = Point.([(171, 171), (185, 111), (202, 109), (202, 189), (328, 160),
(208, 254), (241, 330), (164,252), (69, 278), (139, 208), (72, 148), (168, 172)])
plt = plot(map(a -> a.x, controlpoints), map(a -> a.y, controlpoints)) savefig(plt, "BSplineplot.png")</lang>
Lua
<lang Lua>local function Range(from, to) local range = {} for n = from, to do table.insert(range, n) end return range end
local function Bspline(controlPoints, k) return { controlPoints = controlPoints, n = #controlPoints, k = k, t = Range(1, #controlPoints+k), -- Use a uniform knot vector, delta=1. } end
local function helper(bspline, i, k, x) return (bspline.t[i+k] ~= bspline.t[i]) and (x - bspline.t[i]) / (bspline.t[i+k] - bspline.t[i]) or 0 end
local function calculateBspline(bspline, i, k, x) if k == 1 then return (bspline.t[i] <= x and x < bspline.t[i+1]) and 1 or 0 end return ( helper(bspline, i , k-1, x)) * calculateBspline(bspline, i , k-1, x) + (1-helper(bspline, i+1, k-1, x)) * calculateBspline(bspline, i+1, k-1, x) end
local function round(n) return math.floor(n+.5) end
local function getBsplinePoints(bspline) local points = {}
for x = bspline.t[bspline.k], bspline.t[bspline.n+1]-1 do local sumX = 0 local sumY = 0
for i = 1, bspline.n do local f = calculateBspline(bspline, i, bspline.k, x) sumX = sumX + f * bspline.controlPoints[i].x sumY = sumY + f * bspline.controlPoints[i].y end table.insert(points, {x=round(sumX), y=round(sumY)}) end
return points end
local function Plot(unscaledWidth,unscaledHeight, scaleX,scaleY) local plot = { width = round(unscaledWidth * scaleX), height = round(unscaledHeight * scaleY), scaleX = scaleX, scaleY = scaleY, } for row = 1, plot.height do plot[row] = {} end return plot end
local function raytrace(x,y, x2,y2, visit) local dx = math.abs(x2 - x) local dy = math.abs(y2 - y) local n = 1 + dx + dy local dirX = (x2 > x) and 1 or -1 local dirY = (y2 > y) and 1 or -1 local err = dx - dy dx, dy = 2*dx, 2*dy
for n = 1, n do visit(x, y) if err > 0 then x, err = x+dirX, err-dy else y, err = y+dirY, err+dx end end end
local function plotLine(plot, x1,y1, x2,y2) raytrace(x1,y1, x2,y2, function(x, y) if x >= 0 and y >= 0 and x < plot.width and y < plot.height then plot[y+1][x+1] = true end end) end
local function plotBspline(bspline, plot) if bspline.k > bspline.n or bspline.k < 1 then error("k (= "..bspline.k..") can't be more than "..bspline.n.." or less than 1.") end local points = getBsplinePoints(bspline)
-- Plot the curve. for i = 1, #points-1 do local p1 = points[i] local p2 = points[i+1] plotLine(plot, round(p1.x*plot.scaleX), round(p1.y*plot.scaleY), round(p2.x*plot.scaleX), round(p2.y*plot.scaleY) ) end end
local function printPlot(plot) for row = 1, plot.height do for column = 1, plot.width do io.write(plot[row][column] and "@" or " ") end io.write("\n") end end
local controlPoints = { {x=171, y=171}, {x=185, y=111}, {x=202, y=109}, {x=202, y=189}, {x=328, y=160}, {x=208, y=254}, {x=241, y=330}, {x=164, y=252}, {x= 69, y=278}, {x=139, y=208}, {x= 72, y=148}, {x=168, y=172}, } local k = 4 -- Polynomial degree is one less than this i.e. cubic. local bspline = Bspline(controlPoints, k)
local scaleX = .4 -- Since we print the plot to the console as text let's scale things appropriately. local scaleY = .2 local plot = Plot(350,350, scaleX,scaleY) plotBspline(bspline, plot)
printPlot(plot)</lang>
- Output:
@@@@ @@@@ @@ @@ @@ @@ @@ @@ @@ @@ @@@@@@@@ @@@@@@@@@@@@@@ @@@@@@@@ @@ @@@ @@ @@@ @@ @ @@@ @@ @@ @@ @@@ @ @@ @@ @@@ @@ @@ @@ @@@ @ @@ @@ @@ @@ @@ @@@@@@@ @ @@@@@@@@@@@@@@ @@ @@@@@@@@@ @ @@@@ @ @@@@@ @@ @@@@@ @ @@@@ @@ @@@@@ @ @@@@ @@ @@@
Mathematica/Wolfram Language
<lang Mathematica>Graphics[
BSplineCurve[{{171, 171}, {185, 111}, {202, 109}, {202, 189}, {328, 160}, {208, 254}, {241, 330}, {164, 252}, {69, 278}, {139, 208}, {72, 148}, {168, 172}}, SplineClosed -> True, SplineDegree -> 2]]</lang>
- Output:
Outputs a graphical representation of a B-spline.
Perl
<lang perl>use strict; use warnings; use Class::Struct; use Cairo;
{ package Line;
struct( A => '@', B => '@');
}
my ($WIDTH, $HEIGHT, $W_LINE, $CURVE_F, $DETACHED, $OUTPUT ) =
( 400, 400, 2, 0.25, 0, 'run/b-spline.png' );
my @pt = (
[171, 171], [185, 111], [202, 109], [202, 189], [328, 160], [208, 254], [241, 330], [164, 252], [ 69, 278], [139, 208], [ 72, 148], [168, 172]
); my $cnt = @pt;
sub angle {
my($g) = @_; atan2 $g->B->[1] - $g->A->[1], $g->B->[0] - $g->A->[0]
}
sub control_points {
my($g, $l) = @_;
my $h = Line->new; my $lgt = sqrt( ($g->B->[0] - $l->A->[0])**2 + ($g->B->[1] - $l->A->[1])**2 );
@{$h->B} = @{$l->A}; @{$h->A} = ($g->B->[0] - $lgt * cos(angle $g) , $g->B->[1] - $lgt * sin(angle $g)); my $a = angle $h; my @p1 = ($g->B->[0] + $lgt * cos($a) * $CURVE_F, $g->B->[1] + $lgt * sin($a) * $CURVE_F);
@{$h->A} = @{$g->B}; @{$h->B} = ($l->A->[0] + $lgt * cos(angle $l) , $l->A->[1] + $lgt * sin(angle $l)); $a = angle $h; my @p2 = ($l->A->[0] - $lgt * cos($a) * $CURVE_F, $l->A->[1] - $lgt * sin($a) * $CURVE_F);
\@p1, \@p2
}
my $surf = Cairo::ImageSurface->create ('argb32', $WIDTH, $HEIGHT); my $cr = Cairo::Context->create ($surf); $cr->set_line_width($W_LINE); $cr->move_to($pt[$DETACHED - 1 + $cnt][0], $pt[$DETACHED - 1 + $cnt][1]);
my Line ($g,$l); for my $j ($DETACHED..$cnt-1) {
$g = Line->new( A=>$pt[($j + $cnt - 2) % $cnt], B=>$pt[($j + $cnt - 1) % $cnt]); $l = Line->new( A=>$pt[($j + $cnt + 0) % $cnt], B=>$pt[($j + $cnt + 1) % $cnt]); my($p1,$p2) = control_points($g, $l); $cr->curve_to($$p1[0], $$p1[1], $$p2[0], $$p2[1], $pt[$j][0], $pt[$j][1]);
} $cr->stroke; $surf->write_to_png($OUTPUT);</lang> Output: b-spline.png (offsite image)
Phix
You can run this online here.
-- -- demo\rosetta\B-spline.exw -- ========================= -- -- Use +/- to change the order between k = 1 and k = 4. -- with javascript_semantics include pGUI.e include IupGraph.e constant ctrl_points = {{171, 171}, {185, 111}, {202, 109}, {202, 189}, {328, 160}, {208, 254}, {241, 330}, {164, 252}, { 69, 278}, {139, 208}, { 72, 148}, {168, 172}} integer k = 2, n sequence t function w(integer i, k, x) // B-spline helper function return iff(t[i+k]!=t[i] ? (x-t[i])/(t[i+k]-t[i]) : 0 ) end function function b(integer i, k, x) // B-spline function if k==1 then return iff(t[i]<=x and x<t[i+1] ? 1 : 0) end if return w(i,k-1,x)*b(i,k-1,x) + (1-w(i+1,k-1,x))*b(i+1,k-1,x) end function function b_spline(Ihandle graph) n = length(ctrl_points) t = tagset(n+1+k) // use a uniform knot vector, delta = 1 assert(k<=n+1 and k>=1,"k (= %d) cannot be more than %d or less than 1.",{k,n+1}) sequence px = {}, py = {} for x=t[k] to t[n+1] do atom sumX = 0, sumY = 0 for i=1 to n do atom f = b(i,k,x) sumX += f*ctrl_points[i][1] sumY += f*ctrl_points[i][2] end for px &= round(sumX) py &= round(sumY) end for integer xtick = 40, ytick = 40, xmin = trunc(min(px)/xtick)*xtick, xmax = ceil(max(px)/xtick)*xtick, ymin = trunc(min(py)/ytick)*ytick, ymax = ceil(max(py)/ytick)*ytick IupSetInt(graph,"XTICK",xtick) IupSetInt(graph,"XMIN",xmin) IupSetInt(graph,"XMAX",xmax) IupSetInt(graph,"YTICK",ytick) IupSetInt(graph,"YMIN",ymin) IupSetInt(graph,"YMAX",ymax) sequence graphdata = {{px,py,CD_BLUE}} return graphdata end function procedure set_title(Ihandle dlg) IupSetStrAttribute(dlg, "TITLE", "B-spline curve (order k = %d)",{k}) end procedure function key_cb(Ihandle dlg, atom c) if c=K_ESC then return IUP_CLOSE end if if c='+' then k = min(k+1,4) end if if c='-' then k = max(k-1,1) end if set_title(dlg) IupRedraw(dlg) return IUP_IGNORE end function procedure main() IupOpen() Ihandle graph = IupGraph(b_spline,`RASTERSIZE=600x600`) Ihandle dlg = IupDialog(graph) IupSetCallback(dlg, "KEY_CB", Icallback("key_cb")) set_title(dlg) IupShow(dlg) if platform()!=JS then IupMainLoop() IupClose() end if end procedure main()
Raku
A minimal translation of this C program, by Bernhard R. Fischer. <lang perl6># 20211112 Raku programming solution
use Cairo;
- class point_t { has Num ($.x,$.y) is rw } # get by with two element lists
class line_t { has ($.A,$.B) is rw }
my (\WIDTH, \HEIGHT, \W_LINE, \CURVE_F, \DETACHED, \OUTPUT ) =
400, 400, 2, 0.25, 0, './b-spline.png' ;
my \cnt = #`(Number of points) ( my \pt = [
[171, 171], [185, 111], [202, 109], [202, 189], [328, 160], [208, 254], [241, 330], [164, 252], [ 69, 278], [139, 208], [ 72, 148], [168, 172], ]
).elems;
sub angle(\g) { atan2(g.B.[1] - g.A.[1], g.B.[0] - g.A.[0]) }
sub control_points(\g, \l, @p1, @p2){
- `[ This function calculates the control points. It takes two lines g and l as
* arguments but it takes three lines into account for calculation. This is * line g (P0/P1), line h (P1/P2), and line l (P2/P3). The control points being * calculated are actually those for the middle line h, this is from P1 to P2. * Line g is the predecessor and line l the successor of line h. * @param g Pointer to first line (P0 to P1) * @param l Pointer to third line (P2 to P3) * @param p1 Pointer to memory of first control point. * @param p2 Pointer to memory of second control point. ] my \h = $ = line_t.new;
my \lgt = sqrt([+]([ g.B.[0]-l.A.[0], g.B.[1]-l.A.[1] ]>>²));#length of P1 to P2
h.B = l.A.clone; # end point of 1st tangent # start point of tangent at same distance as end point along 'g' h.A = g.B.[0] - lgt * cos(angle g) , g.B.[1] - lgt * sin(angle g);
my $a = angle h ; # angle of tangent # 1st control point on tangent at distance 'lgt * CURVE_F' @p1 = g.B.[0] + lgt * cos($a) * CURVE_F, g.B.[1] + lgt * sin($a) * CURVE_F;
h.A = g.B.clone; # start point of 2nd tangent # end point of tangent at same distance as start point along 'l' h.B = l.A.[0] + lgt * cos(angle l) , l.A.[1] + lgt * sin(angle l);
$a = angle h; # angle of tangent # 2nd control point on tangent at distance 'lgt * CURVE_F' @p2 = l.A.[0] - lgt * cos($a) * CURVE_F, l.A.[1] - lgt * sin($a) * CURVE_F;
}
given Cairo::Image.create(Cairo::FORMAT_ARGB32, WIDTH, HEIGHT) {
given Cairo::Context.new($_) { my line_t ($g,$l); my (@p1,@p2);
.line_width = W_LINE; .move_to(pt[DETACHED - 1 + cnt].[0], pt[DETACHED - 1 + cnt].[1]);
for DETACHED..^cnt -> \j { $g = line_t.new: A=>pt[(j + cnt - 2) % cnt], B=>pt[(j + cnt - 1) % cnt]; $l = line_t.new: A=>pt[(j + cnt + 0) % cnt], B=>pt[(j + cnt + 1) % cnt]; # Calculate controls points for points pt[j-1] and pt[j]. control_points($g, $l, @p1, @p2); .curve_to(@p1[0], @p1[1], @p2[0], @p2[1], pt[j].[0], pt[j].[1]); } .stroke; }; .write_png(OUTPUT) and die # C return
}</lang>
Output: (Offsite image file)
Wren
In the absence of any clarification on what to use (see Talk page), the following uses a degree of 3 (i.e order k = 4) and a uniform knot vector from 1 to 16 (as there are 12 control points) with a delta of 1.
If one uses a value for k of 1, then the script will simply plot the control points as in the Julia example. <lang ecmascript>import "dome" for Window, Process import "graphics" for Canvas, Color
class BSpline {
construct new(width, height, cpoints, k) { Window.resize(width, height) Canvas.resize(width, height) Window.title = "B-spline curve" _p = cpoints _n = cpoints.count - 1 _k = k _t = (1.._n + 1 + k).toList // use a uniform knot vector, delta = 1 }
// B-spline helper function w(i, k, x) { (_t[i+k] != _t[i]) ? (x - _t[i]) / (_t[i+k] - _t[i]) : 0 } // B-spline function b(i, k, x) { if (k == 1) return (_t[i] <= x && x < _t[i + 1]) ? 1 : 0 return w(i, k-1, x) * b(i, k-1, x) + (1 - w(i+1, k-1, x)) * b(i+1, k-1, x) }
// B-spline points p() { var bpoints = [] for (x in _t[_k-1]..._t[_n + 1]) { var sumX = 0 var sumY = 0 for (i in 0.._n) { var f = b(i, _k, x) sumX = sumX + f * _p[i][0] sumY = sumY + f * _p[i][1] } bpoints.add([sumX.round, sumY.round]) } return bpoints }
init() { if (_k > _n + 1 || _k < 1) { System.print("k (= %(_k)) can't be more than %(_n+1) or less than 1.") Process.exit() } var bpoints = p() // plot the curve for (i in 1...bpoints.count) { Canvas.line(bpoints[i-1][0], bpoints[i-1][1], bpoints[i][0], bpoints[i][1], Color.white) } }
update() {}
draw(alpha) {}
}
var cpoints = [
[171, 171], [185, 111], [202, 109], [202, 189], [328, 160], [208, 254], [241, 330], [164, 252], [ 69, 278], [139, 208], [ 72, 148], [168, 172]
] var k = 4 // polynomial degree is one less than this i.e. cubic var Game = BSpline.new(400, 400, cpoints, k)</lang>