B-spline
B-spline is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
- 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
ALGOL 68
Which is
Suppresses unused parts of the plot.
BEGIN # construct a B-Spline #
# mode to hold a B Spline #
MODE BSPLINE = STRUCT( FLEX[ 1 : 0, 1 : 2 ]INT control points
, INT n
, INT k
, FLEX[ 1 : 0 ]INT t
);
PROC uniform knot vector = ( INT lb, ub )[]INT:
BEGIN
[ lb : ub ]INT range;
FOR n FROM lb TO ub DO range[ n ] := n OD;
range
END # uniform knot vector # ;
PROC calculate bspline = ( REF BSPLINE bs, INT i, k, x )REAL:
IF k = 1
THEN ABS ( ( t OF bs )[ i ] <= x AND x < ( t OF bs )[ i + 1 ] )
ELSE
PROC helper = ( REF BSPLINE bs, INT i, k, x )REAL:
IF ( t OF bs )[ i + k ] /= ( t OF bs )[ i ]
THEN ( x - ( t OF bs )[ i ] ) / ( ( t OF bs )[ i + k ] - ( t OF bs )[ i ] )
ELSE 0
FI # helper # ;
( helper( bs, i, k - 1, x ) ) * calculate bspline( bs, i, k - 1, x )
+ ( 1 - helper( bs, i + 1, k - 1, x ) ) * calculate bspline( bs, i + 1, k - 1, x )
FI # calculate bspline # ;
PROC round = ( REAL n )INT: ENTIER( n + 0.5 );
PROC get bspline points = ( REF BSPLINE bs )[,]INT:
BEGIN
INT p from = ( t OF bs )[ k OF bs ];
INT p to = ( t OF bs )[ 1 + n OF bs ] - 1;
[ p from : p to, 1 : 2 ]INT points;
FOR x FROM p from TO p to DO
REAL sum x := 0;
REAL sum y := 0;
FOR i TO n OF bs DO
REAL f = calculate bspline( bs, i, k OF bs, x );
sum x +:= f * ( control points OF bs )[ i, 1 ];
sum y +:= f * ( control points OF bs )[ i, 2 ]
OD;
points[ x, 1 ] := round( sum x );
points[ x, 2 ] := round( sum y )
OD;
points
END # get bspline points # ;
PROC raytrace = ( INT x0, y0, x2, y2, REF[,]BOOL plot, PROC( REF[,]BOOL, INT, INT )VOID visit )VOID:
BEGIN
INT x := x0;
INT y := y0;
INT dx := ABS ( x2 - x );
INT dy := ABS ( y2 - y );
INT n = 1 + dx + dy;
INT dir x = IF x2 > x THEN 1 ELSE -1 FI;
INT dir y = IF y2 > y THEN 1 ELSE -1 FI;
INT err := dx - dy;
dx *:= 2;
dy *:= 2;
FOR i TO n DO
visit( plot, x, y );
IF err > 0
THEN x +:= dir x; err -:= dy
ELSE y +:= dir y; err +:= dx
FI
OD
END # raytrace # ;
PROC plot line = ( REF[,]BOOL plot, INT x1, y1, x2, y2 )VOID:
raytrace( x1, y1, x2, y2, plot
, ( REF[,]BOOL plot, INT x, INT y )VOID: IF x >= 0
AND y >= 0
AND x < 1 UPB plot
AND y < 2 UPB plot
THEN plot[ x + 1, y + 1 ] := TRUE
FI
);
PROC plot bspline = ( REF BSPLINE bs, REF[,]BOOL plot, REAL scale x, scale y )VOID:
IF k OF bs > n OF bs OR k OF bs < 1 THEN
print( ( "k (= ", whole( k OF bs, 0 ), ") can't be more than ", whole( n OF bs, 0 ) ) );
print( ( " or less than 1.", newline ) );
stop
ELSE
[,]INT points = get bspline points( bs );
# Plot the curve. #
FOR i FROM 1 LWB points TO 1 UPB points - 1 DO
INT p1x = points[ i, 1 ], p1y = points[ i, 2 ];
INT p2x = points[ i + 1, 1 ], p2y = points[ i + 1, 2 ];
plot line( plot
, round( p1x * scale x ), round( p1y * scale y )
, round( p2x * scale x ), round( p2y * scale y )
)
OD
FI # plot bspline # ;
# print the plot - outputs @ or blank depending on whether the point is plotted or not #
PROC print plot = ( [,]BOOL plot )VOID:
FOR row FROM 1 LWB plot
TO BEGIN # find the highest used row #
INT max row := 1 UPB plot;
WHILE IF max row < 1 LWB plot
THEN FALSE
ELSE
BOOL empty row := TRUE;
FOR column FROM 2 LWB plot TO 2 UPB plot
WHILE empty row := NOT plot[ column, max row ]
DO
SKIP
OD;
empty row
FI
DO
max row -:= 1
OD;
max row
END
DO
INT max column := 2 UPB plot;
WHILE IF max column < 2 LWB plot THEN FALSE ELSE NOT plot[ max column, row ] FI
DO
max column -:= 1
OD;
FOR column FROM 2 LWB plot TO max column DO
print( ( IF plot[ column, row ] THEN "@" ELSE " " FI ) )
OD;
print( ( newline ) )
OD # print plot # ;
# task #
[,]INT control 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 )
);
INT k = 4; # Polynomial degree is one less than this i.e. cubic. #
BSPLINE bs := ( control points
, UPB control points
, k
, uniform knot vector( 1, UPB control points + k )
);
REAL scale x = 0.4; # Since we print the plot to the console as text let's scale things appropriately. #
REAL scale y = 0.2;
[ 1 : 350, 1 : 350 ]BOOL plot;
FOR r FROM 1 LWB plot TO 1 UPB plot DO FOR c FROM 2 LWB plot TO 2 UPB plot DO plot[ c, r ] := FALSE OD OD;
plot bspline( bs, plot, scale x, scale y );
print plot( plot )
END- Output:
leading blank lines removed...
@@@@
@@@@
@@
@@
@@
@@
@@
@@
@@
@@
@@@@@@@@
@@@@@@@@@@@@@@
@@@@@@@@
@@
@@@
@@
@@@
@@
@ @@@
@@ @@
@@ @@@
@ @@
@@ @@@
@@ @@
@@ @@@
@ @@
@@ @@
@@ @@
@@@@@@@ @
@@@@@@@@@@@@@@ @@
@@@@@@@@@ @
@@@@ @
@@@@@ @@
@@@@@ @
@@@@ @@
@@@@@ @
@@@@ @@
@@@
Julia
Choose BSpline D of 2, ie degree 1.
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")
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)
- Output:
@@@@
@@@@
@@
@@
@@
@@
@@
@@
@@
@@
@@@@@@@@
@@@@@@@@@@@@@@
@@@@@@@@
@@
@@@
@@
@@@
@@
@ @@@
@@ @@
@@ @@@
@ @@
@@ @@@
@@ @@
@@ @@@
@ @@
@@ @@
@@ @@
@@@@@@@ @
@@@@@@@@@@@@@@ @@
@@@@@@@@@ @
@@@@ @
@@@@@ @@
@@@@@ @
@@@@ @@
@@@@@ @
@@@@ @@
@@@
Mathematica /Wolfram Language
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]]
- Output:
Outputs a graphical representation of a B-spline.
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);
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()
POV-Ray
//cmd: +w640 +h480 +am2 +a0.01
#version 3.7;
global_settings {
assumed_gamma 1
}
//enum CICurve items
#declare _cu_items = array[4]{"Ctrl", "Point", "Curmult", "Matrix"}
#declare _cu_ctrl = 0;
#declare _cu_point = 1;
#declare _cu_curmult = 2;
#declare _cu_matrix = 3;
//enum curve ctrl items
#declare _cu_ctrl_items = array[5]{"DimSeg", "Step", "Len", "Segments", "CurrSeg"}
#declare _cu_dimseg = 0;
#declare _cu_step = 1;
#declare _cu_len = 2;
#declare _cu_segments = 3;
#declare _cu_curseg = 4;
//enum CITCurve items
#declare _cu_tcurve = 0;
#declare _cu_ccurve = 1;
#macro _arrFill(Len, Val)
/*:Create and fill an array of length Len with value Val*/
#local Arr = array[Len];
#for(i,0,Len-1)
#local Arr[i] = Val;
#end
Arr
#end
#macro _CIinit(CIarr, DimSeg, Step, Matrix)
#local Len = dimension_size(CIarr, 1);
#local Lsd = Len - DimSeg;
#if(Lsd < 0)
#error "Too few points in curve array"
#end
#if(mod(Lsd, Step) = 0)
#local Segments = div(Lsd, Step) + 1;
#else
#error "Not the right amount of points in curve array"
#end
#local CurrSeg = -1;
#local CICurve = array[4];
#local CICurve[_cu_ctrl] = array[5]{DimSeg, Step, Len, Segments, CurrSeg};
#local CICurve[_cu_point] = CIarr;
#local CICurve[_cu_curmult] = _arrFill(dimension_size(Matrix,1), <0,0,0>);
#local CICurve[_cu_matrix] = Matrix;
CICurve
#end
#macro _CIdivMatrix(Matrix, Div)
#for(i, 0, dimension_size(Matrix,1)-1)
#for(j, 0, dimension_size(Matrix[i],1)-1)
#local Matrix[i][j] = Matrix[i][j] / Div;
#end
#end
Matrix
#end
#macro CIbspline2(CIarr)
/*:
Quadratic b-spline aka parabolic spline
Three control points define a segment.
The curve does not pass through any control point. It starts half way P0P1
and ends half way P1P2.
The curve is "attracted" to the middel control points.
Number of control points in a multi segment curve is 3 + (n * 1)
*/
#local DimSeg = 3;
#local Step = 1;
#local Div = 2;
#local Matrix = array[3]{
array[3]{ 1,-2, 1},
array[2]{-2, 2},
array[2]{ 1, 1}
}
#local Matrix = _CIdivMatrix(Matrix, Div);
_CIinit(CIarr, DimSeg, Step, Matrix)
#end
#macro CIbspline3(CIarr)
/*: Cubic b-spline
Four control points define a segment.
The curve does not pass through any control point. It goes from near
P1 to near P2
The curve is "attracted" to the middel control points.
Number of control points in a multi segment curve is 4 + (n * 1)
*/
#local DimSeg = 4;
#local Step = 1;
#local Div = 6;
#local Matrix = array[4]{
array[4]{-1, 3,-3, 1},
array[3]{ 3,-6, 3},
array[3]{-3, 0, 3},
array[3]{ 1, 4, 1}
}
#local Matrix = _CIdivMatrix(Matrix, Div);
_CIinit(CIarr, DimSeg, Step, Matrix)
#end
#macro CIget(CICurve, T)
/*: Gets the point on the curve at T
parameter:
CICurve (array, vector): an array of vectors and other parameters as created via one
of the spline initialisation macros.
T (float, [0,1]): position at which to interpolate the curve.
return:
(vector) returns a point on the curve at T.
use example:
#declare TheCurve = CIlinear(PointArray);
#declare P = CIget(TheCurve, 0.1);
*/
/* T goes from 0 to 1 for the whole curve. From current T find the proper curve segment.
Tseg also goes from 0 to 1. Once the proper segment is found, get the proper Tseg.
*/
#local Len = CICurve[_cu_ctrl][_cu_len];
#local Segments = CICurve[_cu_ctrl][_cu_segments];
#local CImult = CICurve[_cu_curmult];
#if(T = 1)
#local Segment = Segments - 1;
#local Tseg = 1;
#else
#local Segment = floor(T * Segments);
#local Tseg = (T * Segments) - Segment;
#end
#if(Segment != CICurve[_cu_ctrl][_cu_curseg])
#local CICurve[_cu_ctrl][_cu_curseg] = Segment;
#local CurrPos = Segment * CICurve[_cu_ctrl][_cu_step];
#local CImatrix = CICurve[_cu_matrix];
#local MaxPos = Len - CICurve[_cu_ctrl][_cu_dimseg] + CICurve[_cu_ctrl][_cu_step];
#if(CurrPos < MaxPos)
#for(i, 0, dimension_size(CImatrix, 1) - 1)
#local CImult[i] = <0, 0, 0>;
#for(j, 0, dimension_size(CImatrix[i], 1) - 1)
#local jSegment = j + CurrPos;
#local CImult[i] = CImult[i] + CImatrix[i][j] * CICurve[_cu_point][jSegment];
#end
#end
#end
#local CICurve[_cu_curmult] = CImult;
#end
#local F = CImult[0];
#for(i, 1, dimension_size(CImult, 1) - 1)
#local F = (F * Tseg) + CImult[i];
#end
F
#end
//== Scene
//spline control points
#declare ControlPoints = array[12]{
<171,171>,
<185,111>,
<202,109>,
<202,189>,
<328,160>,
<208,254>,
<241,330>,
<164,252>,
< 69,278>,
<139,208>,
< 72,148>,
<168,172>
};
//Set up quadratic b-spline (White)
#declare Bspl2 = CIbspline2(ControlPoints)
//get the points on the curve
#for(T, 0, 1, 0.001)
sphere{CIget(Bspl2, T), 2 pigment{rgb 1}}
#end
//Set up cubic b-spline (Blue)
#declare Bspl3 = CIbspline3(ControlPoints)
//get the points on the curve
#for(T, 0, 1, 0.001)
sphere{CIget(Bspl3, T), 2 pigment{rgb <0.2,0.2,1>}}
#end
//control points
sphere{ControlPoints[0] 3 pigment{rgb <1,0,0>}}
#for(i, 1, dimension_size(ControlPoints,1) - 1)
sphere{ControlPoints[i] 3 pigment{rgb <1,0,0>}}
cylinder{ControlPoints[i-1], ControlPoints[i], 0.5 pigment{rgb <1,0,0>}}
#end
camera {
//perspective angle 75
location <190, 220,-250>
right x*image_width/image_height
look_at <190, 220,0>
}
light_source{<0, 0,-3000> color rgb 1}
- Output:
Processing
//Aamrun, 26th June 2022
int abscissae[]={171,185,202,202,328,208,241,164,69,139,72,168};
int ordinates[] = {171,111,109,189 ,160,254 ,330 ,252,278 ,208 ,148 ,172};
void setup() {
size(450, 450);
background(255);
smooth();
noFill();
stroke(0);
beginShape();
for(int i=0;i<abscissae.length;i++){
curveVertex(abscissae[i], ordinates[i]);
}
endShape();
fill(255, 0, 0);
noStroke();
for (int i = 0; i < abscissae.length; i ++) {
ellipse(abscissae[i], ordinates[i], 3, 3);
}
}
Raku
A minimal translation of this C program, by Bernhard R. Fischer.
# 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
}
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.
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)