Ruby script to calc 3D-Spline-Interpolation.

spline_interpolation.rb

#!/usr/local/bin/ruby
#*********************************************
# ３次スプライン補間
# （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
  # ３次スプライン補間
  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