スプライン曲線を描く(Ruby)

複数の点をなめらかにつなぐ「スプライン曲線」を描いてみました。全体的に下のサイトを参考にしました(ありがとうございます!)
Ruby - 3次スプライン補間! - mk-mode BLOG
 
gnuplot で出力してみるとこんな感じです。

 
spline_curve.rb

require 'numo/gnuplot'
require 'matrix'

class SplineCurve
  def initialize(step)
    @step = step    #補完する点の間隔
  end
  
  def generate(points)
    @xs, @ys = points.transpose
    @n = @xs.size - 1
    
    #h, w, v
    h = calc_h
    w = calc_w(h)
    v = calc_v(h, w)
    
    #a, b, c, d
    b = v[0..-2].map {|x| x / 2}
    d = @ys[0..-2]
    a, c = [], []
    @n.times do |i|
      dx = @xs[i + 1] - @xs[i]
      a << (v[i + 1] - v[i]) / (6 * dx)
      c << (@ys[i + 1] - @ys[i]) / dx - dx * (2 * v[i] + v[i + 1]) / 6
    end
    
    Enumerator.new do |yl|
      @xs.each_cons(2).with_index do |xi, i|
        x = xi[0]
        while x < xi[1]
          y = proc do
            dx = x - @xs[i]
            a[i] * dx ** 3 + b[i] * dx ** 2 + c[i] * dx + d[i]
          end
          
          yl << [x, y.call]
          x += @step
        end
        x = xi[1]
        yl << [x, y.call]
      end
    end
  end
  
  def calc_h
    (0...@n).map {|i| @xs[i + 1] - @xs[i]}
  end
  
  def calc_w(h)
    (1...@n).map do |i|
      6 * ((@ys[i + 1] - @ys[i]) / h[i] - (@ys[i] - @ys[i - 1]) / h[i - 1])
    end
  end
  
  def calc_v(h, w)
    m = Array.new(@n - 1) {[0] * (@n - 1)}
    (@n - 1).times do |i|
      m[i][i] = 2 * (h[i] + h[i + 1])
      next if i.zero?
      m[i - 1][i] = h[i]
      m[i][i - 1] = h[i]
    end
    [0] + Matrix[*m].lup.solve(w).to_a + [0]
  end
end


if __FILE__ == $0
  points = [[0.0, 0.8], [2.0, 3.2], [3.0, 2.8], [5.0, 4.5], [7.0, 2.5], [8.0, 1.9]]
  x, y = SplineCurve.new(0.1).generate(points).to_a.transpose
  
  Numo.gnuplot do
    set title: "Spline Curve"
    set terminal: "png"
    set output: "spline_curve.png"
    unset :key
    plot [x, y, w: :l],
         [*points.transpose, w: :points, pt: 7, lc: [rgb: "orange"], pointsize: 2]
  end
end

points が点の座標を入れた配列を与えます。SplineCurve#generate はスプライン補完された点の座標 [x, y] を順に与える Enumerator を返します。ここでは to_a して一気に配列にしているので、あまり意味はないですが。