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 cbs, INT cbi, cbk, cbx )REAL:
IF cbk = 1
THEN ABS ( ( t OF cbs )[ cbi ] <= cbx AND cbx < ( t OF cbs )[ cbi + 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( cbs, cbi, cbk - 1, cbx ) ) * calculate bspline( cbs, cbi, cbk - 1, cbx )
+ ( 1 - helper( cbs, cbi + 1, cbk - 1, cbx ) ) * calculate bspline( cbs, cbi + 1, cbk - 1, cbx )
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 rtplot, 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;
TO n DO
visit( rtplot, 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 plplot, INT x1, y1, x2, y2 )VOID:
raytrace( x1, y1, x2, y2, plplot
, ( 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 bsplot, 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( bsplot
, 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 # ;
BEGIN # 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
END- Output:
leading blank lines removed...
@@@@
@@@@
@@
@@
@@
@@
@@
@@
@@
@@
@@@@@@@@
@@@@@@@@@@@@@@
@@@@@@@@
@@
@@@
@@
@@@
@@
@ @@@
@@ @@
@@ @@@
@ @@
@@ @@@
@@ @@
@@ @@@
@ @@
@@ @@
@@ @@
@@@@@@@ @
@@@@@@@@@@@@@@ @@
@@@@@@@@@ @
@@@@ @
@@@@@ @@
@@@@@ @
@@@@ @@
@@@@@ @
@@@@ @@
@@@
FreeBASIC
Type Point2D
x As Integer
y As Integer
End Type
Type BSpline
controlPoints(Any) As Point2D
n As Integer
k As Integer
t(Any) As Integer
End Type
Function CreateRange(fromVal As Integer, toVal As Integer) As Integer Ptr
Dim As Integer Ptr range = CAllocate((toVal - fromVal + 1) * SizeOf(Integer))
For i As Integer = 0 To toVal - fromVal
range[i] = fromVal + i
Next
Return range
End Function
Function CreateBspline(controlPoints() As Point2D, k As Integer) As BSpline
Dim bs As BSpline
nn = UBound(controlPoints)
kk = k
ReDim bs.controlPoints(nn)
For i As Integer = 0 To nn
bs.controlPoints(i) = controlPoints(i)
Next
ReDim bs.t(nn + k)
For i As Integer = 0 To nn + k
bs.t(i) = i + 1
Next
Return bs
End Function
Function Helper(bs As BSpline, i As Integer, k As Integer, x As Single) As Single
If bs.t(i + k) <> bs.t(i) Then Return (x - bs.t(i)) / (bs.t(i + k) - bs.t(i))
Return 0
End Function
Function CalculateBspline(bs As BSpline, i As Integer, k As Integer, x As Single) As Single
If k = 1 Then
If bs.t(i) <= x AndAlso x < bs.t(i + 1) Then Return 1
Return 0
End If
Return Helper(bs, i, k-1, x) * CalculateBspline(bs, i, k-1, x) + _
(1 - Helper(bs, i+1, k-1, x)) * CalculateBspline(bs, i+1, k-1, x)
End Function
Sub PlotBspline(bs As BSpline)
If kk > nn OrElse kk < 1 Then
Print "k (= " & kk & ") can't be more than " & nn & " or less than 1."
Exit Sub
End If
Dim As Integer lastX, lastY
Dim As Boolean firstPoint = True
ScreenRes 400, 400, 32
WindowTitle "B-spline curve"
Cls
' Draw control points
For i As Integer = 0 To nn
Circle (bs.controlPoints(i).x, bs.controlPoints(i).y), 3, RGB(255,0,0)
Next
' Calculate and draw B-spline points
For x As Single = bs.t(kk-1) To bs.t(nn + 1) Step 0.1
Dim As Single sumX = 0, sumY = 0
For i As Integer = 0 To nn
Dim As Single f = CalculateBspline(bs, i, kk, x)
sumX += f * bs.controlPoints(i).x
sumY += f * bs.controlPoints(i).y
Next
If firstPoint Then
lastX = CInt(sumX)
lastY = CInt(sumY)
firstPoint = False
Else
Line (lastX, lastY)-(CInt(sumX), CInt(sumY)), RGB(255,255,255)
lastX = CInt(sumX)
lastY = CInt(sumY)
End If
Next
End Sub
' Main program
Dim controlPoints(11) As Point2D
controlPoints(0) = Type<Point2D>(171, 171)
controlPoints(1) = Type<Point2D>(185, 111)
controlPoints(2) = Type<Point2D>(202, 109)
controlPoints(3) = Type<Point2D>(202, 189)
controlPoints(4) = Type<Point2D>(328, 160)
controlPoints(5) = Type<Point2D>(208, 254)
controlPoints(6) = Type<Point2D>(241, 330)
controlPoints(7) = Type<Point2D>(164, 252)
controlPoints(8) = Type<Point2D>(69, 278)
controlPoints(9) = Type<Point2D>(139, 208)
controlPoints(10) = Type<Point2D>(72, 148)
controlPoints(11) = Type<Point2D>(168, 172)
Dim bs As BSpline = CreateBspline(controlPoints(), 4)
PlotBspline(bs)
Bsave "B-spline_FreeBASIC.bmp", 0
Sleep
Java
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Paths;
import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;
public final class BSpline {
public static void main(String[] args) throws IOException {
List<Point> controlPoints = List.of(
new Point(171.0, 171.0), new Point(185.0, 111.0), new Point(202.0, 109.0),
new Point(202.0, 189.0), new Point(328.0, 160.0), new Point(208.0, 254.0),
new Point(241.0, 330.0), new Point(164.0, 252.0), new Point( 69.0, 278.0),
new Point(139.0, 208.0), new Point( 72.0, 148.0), new Point(168.0, 172.0) );
List<Double> steps = IntStream.range(0, 1_001).mapToObj( i -> (double) i * 0.001 ).toList();
int degree = 2; // The degree of the spline curve
List<Double> knots = IntStream.range(0, controlPoints.size() + degree + 1)
.mapToObj(Double::valueOf).toList();
List<Point> bSpline2 = bSpline(steps, degree, controlPoints, knots);
degree = 3;
knots = IntStream.range(0, controlPoints.size() + degree + 1)
.mapToObj(Double::valueOf).toList();
List<Point> bSpline3 = bSpline(steps, degree, controlPoints, knots);
final int svgSize = 350;
StringBuilder text = new StringBuilder();
text.append("<svg xmlns='http://www.w3.org/2000/svg'");
text.append(" width='" + svgSize + "' height='" + svgSize + "'>\n");
text.append("<rect width='100%' height='100%' fill='black'/>\n");
for ( int i = 0; i < controlPoints.size() - 1; i++ ) {
text.append( // Join controlPoints with lines in red
"<line x1='" + controlPoints.get(i).x + "' y1='" + ( svgSize - controlPoints.get(i).y )
+ "' x2='" + controlPoints.get(i + 1).x + "' y2='" + ( svgSize - controlPoints.get(i + 1).y )
+ "' stroke='red' stroke-width='1'/>");
}
for ( int i = 0; i < bSpline2.size() - 1; i++ ) {
text.append( // Draw bSpline degree 2 curve in white
"<line x1='" + bSpline2.get(i).x + "' y1='" + ( svgSize - bSpline2.get(i).y )
+ "' x2='" + bSpline2.get(i + 1).x + "' y2='" + ( svgSize - bSpline2.get(i + 1).y )
+ "' stroke='white' stroke-width='3'/>");
}
for ( int i = 0; i < bSpline3.size() - 1; i++ ) {
text.append( // Draw bSpline degree 3 curve in cyan
"<line x1='" + bSpline3.get(i).x + "' y1='" + ( svgSize - bSpline3.get(i).y )
+ "' x2='" + bSpline3.get(i + 1).x + "' y2='" + ( svgSize - bSpline3.get(i + 1).y )
+ "' stroke='cyan' stroke-width='3'/>");
}
text.append("'/>\n</svg>\n");
Files.write(Paths.get("./bSplineJava.svg"), text.toString().getBytes());
}
/**
* Return a list of points belonging to the bSpline curve formed from the given control points.
* Implements de Boor's algorithm; for information visit
* https://en.wikipedia.org/wiki/De_Boor%27s_algorithm
*/
public static List<Point> bSpline(
List<Double> steps, int degree, List<Point> controlPoints, List<Double> knots) {
if ( degree < 1 || degree >= controlPoints.size() ) {
throw new AssertionError("Degree must be > 1 and < controlPoints.size(): " + degree);
}
if ( knots.size() != controlPoints.size() + degree + 1 ) {
throw new AssertionError("Invalid knot vector length: " + knots.size());
}
List<Point> bSplinePoints = new ArrayList<Point>(steps.size());
final double low = knots.get(degree);
final double high = knots.get(controlPoints.size());
steps.forEach( step -> {
// Transform the step values so that they match the range of values of the spline sections
step = step * ( high - low ) + low;
if ( step < low || step > high ) {
throw new AssertionError("Invalid step value: " + step);
}
// Locate the spline segment corresponding to the current step value
int splineSegment = degree;
while ( splineSegment < controlPoints.size() &&
( step < knots.get(splineSegment) || step > knots.get(splineSegment + 1) ) ) {
splineSegment += 1;
}
// Copy 'controlPoints' on each iteration because they are mutated by the algorithm
List<Point> copyPoints = new ArrayList<Point>(controlPoints);
// de Boor's algorithm
for ( int level = 1; level <= degree + 1; level++ ) {
for ( int i = splineSegment; i > splineSegment - degree - 1 + level; i-- ) {
final double alpha = ( step - knots.get(i) ) /
( knots.get(i + degree + 1 - level) - knots.get(i) );
final double x = ( 1.0 - alpha ) * copyPoints.get(i - 1).x + alpha * copyPoints.get(i).x;
final double y = ( 1.0 - alpha ) * copyPoints.get(i - 1).y + alpha * copyPoints.get(i).y;
copyPoints.set(i, new Point(x, y));
}
}
bSplinePoints.addLast(copyPoints.get(splineSegment));
} );
return bSplinePoints;
}
private static record Point(double x, double y) {}
}
- Output:
JavaScript
// ===== Points supplied by the task. =====
const points = [
[171, 171],
[185, 111],
[202, 109],
[328, 160],
[208, 254],
[241, 330],
[164, 252],
[69, 278],
[139, 208],
[72, 148],
[168, 172]
];
// ===== Helper functions =====
const lerp = (a, b, t) => a + (b - a) * t;
const lerpPoint = (a, b, t) => a.map((a, i) => lerp(a, b[i], t));
const computeMidpoint = (a, b) => lerpPoint(a, b, 1 / 2);
const computeOneThird = (a, b) => lerpPoint(a, b, 1 / 3);
const computeTwoThird = (a, b) => lerpPoint(a, b, 2 / 3);
// ===== Create a canvas to draw the spline. =====
const canvas = document.createElement("canvas");
// Since no image size was given, it is assumed to be 500 x 500.
canvas.width = 500;
canvas.height = 500;
const ctx = canvas.getContext("2d");
if (!ctx) throw new Error("Canvas 2D API is not supported on this device.");
const [[firstX, firstY], ...nextPoints] = points;
// ===== Draws a 2px black polyline connecting all control points. =====
ctx.beginPath();
ctx.moveTo(firstX, firstY);
for (const [x, y] of nextPoints) ctx.lineTo(x, y);
ctx.lineWidth = 2;
ctx.strokeStyle = "black";
ctx.stroke();
ctx.closePath();
// ===== Draws spline specified in the task. =====
const [first, second] = points;
const [last, secondLast] = points.slice(-2).reverse();
ctx.beginPath();
drawStartingCurve(first, second);
for (const i of Array(points.length - 2).keys()) {
const [a, b, c] = points.slice(i, i + 3);
drawMiddleCurve(a, b, c);
}
drawEndingCurve(last, secondLast);
ctx.lineWidth = 2;
ctx.strokeStyle = "red";
ctx.stroke();
ctx.closePath();
// ==== Creates a link to allow the canvas image to be downloaded. ====
const dataURL = canvas.toDataURL("image/png");
const link = document.createElement("a");
link.href = dataURL;
link.download = "b-spline-js.png";
link.click();
// ==== Spline drawing functions =====
function drawStartingCurve(first, second) {
const oneThird = computeOneThird(first, second);
const [firstX, firstY] = first;
const [oneThirdX, oneThirdY] = oneThird;
const [oneSixthX, oneSixthY] = computeMidpoint(first, oneThird)
ctx.moveTo(firstX, firstY);
ctx.quadraticCurveTo(oneThirdX, oneThirdY, oneSixthX, oneSixthY);
}
function drawMiddleCurve(a, b, c) {
const oneThirdAB = computeOneThird(a, b);
const twoThirdAB = computeTwoThird(a, b);
const oneThirdBC = computeOneThird(b, c);
const twoThirdAB_OneThirdBC_Midpoint = computeMidpoint(twoThirdAB, oneThirdBC);
const [oneThirdAB_X, oneThirdAB_Y] = oneThirdAB;
const [twoThirdAB_X, twoThirdAB_Y] = twoThirdAB;
const [twoThirdAB_OneThirdBC_MidpointX, twoThirdAB_OneThirdBC_MidpointY] = twoThirdAB_OneThirdBC_Midpoint;
ctx.bezierCurveTo(oneThirdAB_X, oneThirdAB_Y, twoThirdAB_X, twoThirdAB_Y, twoThirdAB_OneThirdBC_MidpointX, twoThirdAB_OneThirdBC_MidpointY);
}
function drawEndingCurve(last, secondLast) {
const oneThird = computeOneThird(secondLast, last);
const twoThird = computeTwoThird(secondLast, last);
const [lastX, lastY] = last;
const [oneThirdX, oneThirdY] = oneThird;
const [fiveSixthsX, fiveSixthsY] = computeMidpoint(twoThird, last);
ctx.quadraticCurveTo(oneThirdX, oneThirdY, fiveSixthsX, fiveSixthsY);
ctx.lineTo(lastX, lastY);
}
- Output:
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.
Nim
import std/[math, strformat]
import cairo
type
Point = tuple[x, y: float]
BSpline = object
controlPoints: seq[Point]
n, k: int
t: seq[float]
proc initBSpline(controlPoints: openArray[Point]; k: int): BSpline =
result.controlPoints = @controlPoints
result.n = controlPoints.high
result.k = k
for i in 0..(result.n + k):
result.t.add (i + 1).toFloat
proc helper(bs: BSpline; i, k: int; x: float): float =
if bs.t[i + k] != bs.t[i]:
return (x - bs.t[i]) / (bs.t[i + k] - bs.t[i])
proc calculate(bs: BSpline; i, k: int; x: float): float =
result = if k == 1:
if bs.t[i] <= x and x < bs.t[i + 1]: 1
else: 0
else:
bs.helper(i, k - 1, x) * bs.calculate(i, k - 1, x) +
(1 - bs.helper(i + 1, k - 1, x)) * bs.calculate(i + 1, k - 1, x)
proc plot(bs: BSpline) =
if bs.k > bs.n or bs.k < 1:
quit &"k (= {bs.k}) can't be more than {bs.n} or less than 1.", QuitFailure
var lastx, lasty: float
var firstPoint = true
let surface = imageSurfaceCreate(FormatRgb24, 400, 400)
let context = surface.create()
context.setSourceRgb(0.0, 0.0, 0.0)
context.paint()
# Draw control points.
context.setSourceRgb(1.0, 0.0, 0.0)
for cp in bs.controlPoints:
context.arc(cp.x, cp.y, 3, 0, TAU)
context.stroke()
# Calculate and draw B-spline points.
context.setSourceRgb(1.0, 1.0, 1.0)
var x = bs.t[bs.k - 1]
while x <= bs.t[bs.n]:
var sumx, sumy = 0.0
for i, cp in bs.controlPoints:
let f = bs.calculate(i, bs.k, x)
sumx += f * cp.x
sumy += f * cp.y
if firstPoint:
lastx = sumx
lasty = sumy
firstPoint = false
else:
context.moveTo(lastx, lasty)
context.lineTo(sumx, sumy)
context.stroke()
lastx = sumx
lasty = sumy
x += 0.1
if surface.writeToPng("bspline.png") != StatusSuccess:
quit "Error while saving file", QuitFailure
surface.destroy()
const ControlPoints = [Point (171.0, 171.0), (185.0, 111.0),
(202.0, 109.0), (202.0, 189.0),
(328.0, 160.0), (208.0, 254.0),
(241.0, 330.0), (164.0, 252.0),
( 69.0, 278.0), (139.0, 208.0),
( 72.0, 148.0), (168.0, 172.0)]
let bspline = initBSpline(ControlPoints, 4)
bspline.plot()
- Output:
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);
}
}
QB64
Type Point2D
x As Integer
y As Integer
End Type
Screen 12
Window (0, 0)-(400, 400)
Dim controlPoints(0 To 100) As Point2D
Dim t(0 To 100) As Integer
Dim k As Integer: k = 4
Dim n As Integer: n = 11
' Initialize control points
Data 171,171,185,111,202,109,202,189,328,160,208,254
Data 241,330,164,252,69,278,139,208,72,148,168,172
For i = 0 To n: Read controlPoints(i).x, controlPoints(i).y: Next i
For i = 0 To n + k: t(i) = i + 1: Next i
PlotBspline n, k, t(), controlPoints()
While InKey$ = "": Wend
End
Function CalculateBspline! (t() As Integer, i As Integer, k As Integer, x As Single)
If k = 1 Then
If t(i) <= x And x < t(i + 1) Then
CalculateBspline! = 1
Else
CalculateBspline! = 0
End If
Exit Function
End If
CalculateBspline! = Helper!(t(), i, k - 1, x) * CalculateBspline!(t(), i, k - 1, x) + (1 - Helper!(t(), i + 1, k - 1, x)) * CalculateBspline!(t(), i + 1, k - 1, x)
End Function
Function Helper! (t() As Integer, i As Integer, k As Integer, x As Single)
If t(i + k) <> t(i) Then
Helper! = (x - t(i)) / (t(i + k) - t(i))
Else
Helper! = 0
End If
End Function
Sub PlotBspline (n As Integer, k As Integer, t() As Integer, controlPoints() As Point2D)
If k > n Or k < 1 Then
Print "k (="; k; ") can't be more than"; n; "or less than 1."
Exit Sub
End If
Cls
' Draw control points
For i = 0 To n: Circle (controlPoints(i).x, controlPoints(i).y), 3, 4: Next i
Dim lastX As Integer, lastY As Integer
Dim firstPoint As Integer: firstPoint = 1
' Calculate and draw B-spline points
x! = t(k - 1)
Do While x! <= t(n + 1)
sumX! = 0
sumY! = 0
For i = 0 To n
f! = CalculateBspline!(t(), i, k, x!)
sumX! = sumX! + f! * controlPoints(i).x
sumY! = sumY! + f! * controlPoints(i).y
Next i
If firstPoint = 1 Then
lastX = Int(sumX! + .5)
lastY = Int(sumY! + .5)
firstPoint = 0
Else
Line (lastX, lastY)-(Int(sumX! + .5), Int(sumY! + .5)), 15
lastX = Int(sumX! + .5)
lastY = Int(sumY! + .5)
End If
x! = x! + .1
Loop
End Sub
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)