Skip to content

Instantly share code, notes, and snippets.

@komasaru
Last active April 19, 2018 06:01
Show Gist options
  • Save komasaru/01b9a98af2a6cae0a5df to your computer and use it in GitHub Desktop.
Save komasaru/01b9a98af2a6cae0a5df to your computer and use it in GitHub Desktop.
Ruby script to calc 3D-Spline-Interpolation.
#!/usr/local/bin/ruby
#*********************************************
# 3次スプライン補間
# (gnuplot によるグラフ描画付き)
#*********************************************
#
require 'gnuplot'
# N + 1 個の点
X = [0.0, 2.0, 3.0, 5.0, 7.0, 8.0]
Y = [0.8, 3.2, 2.8, 4.5, 2.5, 1.9]
class SplineInterpolation
def initialize
@x, @y = X, Y
@n = @x.size - 1
h = calc_h
w = calc_w(h)
matrix = gen_matrix(h, w)
v = [0, gauss_jordan(matrix), 0].flatten
@b = calc_b(v)
@a = calc_a(v)
@d = calc_d(v)
@c = calc_c(v)
end
def interpolate(t)
i = search_i(t)
return @a[i] * (t - @x[i]) ** 3 \
+ @b[i] * (t - @x[i]) ** 2 \
+ @c[i] * (t - @x[i]) \
+ @d[i]
rescue => e
$stderr.puts "[#{e.class}] #{e.message}"
exit 1
end
private
def calc_h
return (0..@n - 1).map { |i| @x[i + 1] - @x[i] }
end
def calc_w(h)
return (1..@n -1).map do |i|
6 * ((@y[i + 1] - @y[i]) / h[i] - (@y[i] - @y[i - 1]) / h[i - 1])
end
end
def gen_matrix(h, w)
matrix = Array.new(@n - 1).map { Array.new(@n, 0) }
0.upto(@n - 2) do |i|
matrix[i][i] = 2 * (h[i] + h[i + 1])
matrix[i][-1] = w[i]
next if i == 0
matrix[i - 1][i] = h[i]
matrix[i][i - 1] = h[i]
end
return matrix
end
def gauss_jordan(matrix)
v = Array.new
n = @n - 1
0.upto(n - 1) do |k|
p = matrix[k][k]
k.upto(n) { |j| matrix[k][j] /= p.to_f }
0.upto(n - 1) do |i|
next if i == k
d = matrix[i][k]
k.upto(n) { |j| matrix[i][j] -= d * matrix[k][j] }
end
end
matrix.each { |row| v << row[-1] }
return v
end
def calc_a(v)
return (0..@n - 1).map { |i| (v[i + 1] - v[i]) / (6 * (@x[i + 1] - @x[i])) }
end
def calc_b(v)
return (0..@n - 1).map { |i| v[i] / 2.0 }
end
def calc_c(v)
return (0..@n - 1).map do |i|
(@y[i + 1] - @y[i]) / (@x[i + 1] - @x[i]) \
- (@x[i + 1] - @x[i]) * (2 * v[i] + v[i + 1]) / 6
end
end
def calc_d(v)
return @y
end
def search_i(t)
i, j = 0, @x.size - 1
while i < j
k = (i + j) / 2
if @x[k] < t
i = k + 1
else
j = k
end
end
i -= 1 if i > 0
return i
end
end
class Graph
def initialize(x_g, y_g)
@x_g, @y_g, @x, @y = x_g, y_g, X, Y
end
def plot
Gnuplot.open do |gp|
Gnuplot::Plot.new(gp) do |plot|
plot.terminal "png enhanced font 'IPA P ゴシック' fontscale 1.2"
plot.output "spline_interpolation.png"
plot.title "スプライン補間"
plot.xlabel "x"
plot.ylabel "y"
plot.grid
# 計算によって得られた点
plot.data << Gnuplot::DataSet.new([@x_g, @y_g]) do |ds|
ds.with = "points"
ds.linewidth = 2
ds.linecolor = 3
ds.notitle
end
# 予め与えらた点
plot.data << Gnuplot::DataSet.new([@x, @y]) do |ds|
ds.with = "points"
ds.linewidth = 2
ds.linecolor = 1
ds.notitle
end
end
end
rescue => e
$stderr.puts "[#{e.class}] #{e.message}"
exit 1
end
end
if __FILE__ == $0
# グラフ用配列
x_g, y_g = Array.new, Array.new
# 3次スプライン補間
obj = SplineInterpolation.new
X[0].step(X[-1], 0.1) do |x|
y = obj.interpolate(x)
puts "%8.4f, %8.4f" % [x, y]
x_g << x
y_g << y
end
# グラフ描画
Graph.new(x_g, y_g).plot
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment