B-spline: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎{{header|ALGOL 68}}: Use spaces instead of tabs, tweak)
No edit summary
 
Line 566: Line 566:
<span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
<!--</syntaxhighlight>-->


=={{header|POV-Ray}}==

<syntaxhighlight lang="pov">

//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}


</syntaxhighlight>
{{out}}
[[File:BsplinePOVRay.jpg]]


=={{header|Processing}}==
=={{header|Processing}}==
<syntaxhighlight lang="java">
<syntaxhighlight lang="java">

Latest revision as of 07:14, 3 July 2024

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

Translation of: Lua

Which is

Translation of: Wren

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

Translation of: Wren
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

Translation of: Raku
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

Translation of: Wren
Library: Phix/pGUI
Library: Phix/online

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

Library: DOME

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)
Retrieved from "https://rosettacode.org/wiki/B-spline?oldid=366020"