This is a description of a class of cubic spline functions that have the
property that the endpoints of the curve will pass through the first and
last control points, and will pass
near the interior control points.
I've provided an implementation in ruby (see bottom of this page).
Let {Pi : i=1,N; N>= 4} be the set of control points
Let P0 = P1 and PN+1 = PN
Draw the curve for each Segment Si, i=1, N-1, using
Si(u) = Σ[j=1,N]Bij(u)Pi-2+j
0 <= u < 1
The B
ij are
blending functions defined as:
First segment:
B11 = (1-u)3
B12 = (7u3 - 18u2 + 12u) / 4
B13 = (-11u3 + 18u2) / 12
B14 = u3 / 6
Second segment:
B21 = (1 - u)3 / 4
B22 = (7u3 - 15u2 + 3u + 7) / 12
B23 = (3u3 - 3u2 + 3u + 1) / 6
B24 = u3 / 6
Middle segments:
B31 = (1 - u)3 / 6
B32 = (3u3 - 6u2 + 4) / 6
B33 = (3(-u3 + u2 + u) + 1) / 6
B34 = u3 / 6
Next-to-Last segment:
B(N-1)j = B2(5-j)(1 - u)
Last segment:
BNj = B1(5-j)(1 - u)
In practice you select a set of u
k and draw straight lines between them.
By selecting the number of such points to be large enough (10 is more than enough for most
applications) the curve will look good. This also says that you can pre-compute
the B
ij(u
k) and cache those values in a table, for better performance (it matters, believe me).
If you're wondering where the magic-number coefficients came from: they are selected to ensure that the line segments match at their end-points and that their first-derivatives also match (i.e. if we treat u as continuous rather than a discrete set of u
k values, the derivatives would match). The following table shows the values for u=0 and 1 for each of the 20 B
ij values.
| 0 | 1 |
|---|
| *B1j * | - |
|---|
| 1 | 0 |
| 0 | 3/12 |
| 0 | 7/12 |
| 0 | 2/12 |
| *B2j * | - |
|---|
| 3/12 | 0 |
| 7/12 | 2/12 |
| 2/12 | 8/12 |
| 0 | 2/12 |
| *B3j * | - |
|---|
| 2/12 | 0 |
| 8/12 | 2/12 |
| 2/12 | 8/12 |
| 0 | 2/12 |
| *B4j * | - |
|---|
| 2/12 | 0 |
| 8/12 | 2/12 |
| 2/12 | 7/12 |
| 0 | 2/12 |
| *B5j * | - |
|---|
| 2/12 | 0 |
| 7/12 | 0 |
| 3/12 | 0 |
| 0 | 1 |
class Point
attr_accessor :x, :y
def initialize(x, y)
@x, @y = x, y
end
def to_s
"(#{@x},#{@y})"
end
end
class BSpline
attr_reader :uvalues, :b
attr_accessor :points
def initialize(points, uvalues = 10)
@points = points
if uvalues.kind_of? Array
@uvalues = uvalues
elsif uvalues.kind_of? Numeric
@uvalues = []
d = uvalues.to_f
0.upto(uvalues) {|i| @uvalues << i / d}
end
recalc_bii
end
# define the basic Bxx used for defining the final Bij.
# Each Bxx is a set of coefficents as follows:
# [u^0,u^1,u^2,u^3,(1-u)^3,divisor]
# where the divisor is applied to the whole expression.
# The coefficients can be entered as integers but will be
# evaluated as floating point.
BXY = {:b00 => [0, 0, 0, 0, 1, 1],
:b01 => [0, 12,-18, 7, 0, 4],
:b02 => [0, 0, 18,-11, 0, 12],
:b03 => [0, 0, 0, 1, 0, 6],
:b10 => [0, 0, 0, 0, 1, 4],
:b11 => [7, 3,-15, 7, 0, 12],
:b12 => [1, 3, -3, 3, 0, 6],
:b13 => [0, 0, 0, 1, 0, 6],
:b20 => [0, 0, 0, 0, 1, 6],
:b21 => [4, 0, -6, 3, 0, 6],
:b22 => [1, 3, 3, -3, 0, 6],
:b23 => [0, 0, 0, 1, 0, 6]}
def eval_points
ret = []
0.upto(@uvalues.length - 2) {|k| ret << eval_point([0, 0, 1, 2], 0, k) }
0.upto(@uvalues.length - 2) {|k| ret << eval_point([0, 1, 2, 3], 1, k) }
np = @points.length
1.upto(np - 5) do |j|
0.upto(@uvalues.length - 2) {|k| ret << eval_point([j, j+1, j+2, j+3], 2, k)}
end
if np > 4
0.upto(@uvalues.length - 2) {|k| ret << eval_point([np-4, np-3, np-2, np-1], 3, k)}
end
0.upto(@uvalues.length - 1) {|k| ret << eval_point([np-3,np-2,np-1,np-1], 4, k)}
ret
end
private
def eval_point(indices, i_base, k)
x, y = 0, 0
0.upto(3) do |j|
x += @points[indices[j]].x * @b[i_base][j][k]
y += @points[indices[j]].y * @b[i_base][j][k]
end
Point.new(x, y)
end
def eval_base_bxy(a_sym, u1, u2, u3, uinvcube)
coeff = BXY[a_sym].map{|z| z.to_f}
(coeff[0] + coeff[1] * u1 + coeff[2] * u2 + coeff[3] * u3 + coeff[4] * uinvcube) / coeff[5]
end
def recalc_bii
z = [-1] * @uvalues.length
@b = []
0.upto(4) do |i|
y = []
0.upto(3) {|j| y << z.clone}
@b << y.clone
end
@uvalues.each_with_index do |u1, k|
u2 = u1 * u1
u3 = u2 * u1
uinv = 1.0 - u1
uinvcube = uinv * uinv * uinv
0.upto 2 do |i|
0.upto 3 do |j|
@b[i][j][k] = eval_base_bxy(['b', i.to_s, j.to_s].join('').to_sym, u1, u2, u3, uinvcube)
end
end
u1, uinv = uinv, u1
u2 = u1 * u1
u3 = u2 * u1
uinvcube = uinv * uinv * uinv
3.upto(4) do |i|
0.upto 3 do |j|
@b[i][j][k] = eval_base_bxy(['b', (4-i).to_s, (3-j).to_s].join('').to_sym, u1, u2, u3, uinvcube)
end
end
end
end
end
if __FILE__ == $0
b = BSpline.new([Point.new(0, 0), Point.new(0, 2), Point.new(-1, 1),
Point.new(0, 3), Point.new(1,1), Point.new(5,0)], 5)
puts b.points.join(";")
puts b.eval_points.join("\n")
end
Books RSS
