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

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

 
spline_curve.rb

require 'numo/gnuplot'
require 'matrix'

class SplineCurve
  def initialize
    @step = 0.1
  end
  attr_writer :step
  
  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|
      a << (v[i + 1] - v[i]) / (6 * (@xs[i + 1] - @xs[i]))
      c << ((@ys[i + 1] - @ys[i]) / (@xs[i + 1] - @xs[i]) -
              (@xs[i + 1] - @xs[i]) * (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
            a[i] * (x - @xs[i]) ** 3 + b[i] * (x - @xs[i]) ** 2 +
                c[i] * (x - @xs[i]) + d[i]
          end
          
          yl << [x, y.call]
          x += @step
        end
        x = xi[1]
        yl << [x, y.call]
      end
    end
  end
  
  def calc_h
    0.upto(@n - 1).with_object([]) {|i, h| h << (@xs[i + 1] - @xs[i])}
  end
  
  def calc_w(h)
    1.upto(@n - 1).with_object([]) do |i, w|
      w << 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.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 して一気に配列にしているので、あまり意味はないですが。

なお、@step は補完する点の間隔(dx)を表しています。setter(attr_writer 指定してあります)で変更することができます。