#!ruby -Ks
# -*- mode: ruby; coding: sjis -*-
# Last updated: <2016/02/14 05:39:09 +0900>
#
# DXRubyで3次スプライン曲線を描画
# xが右方向に増えていくことが前提、らしい？
#
# 参考ページ
#
# [PHP]3次スプライン曲線を使ったスプライン補間 | PHP Archive
# http://php-archive.net/php/cubic-spline/
#
# 数値計算
# http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/num/num.htm#7.2

require 'dxruby'

# 3次スプライン補間
class CubicSpline

  # コンストラクタ
  # @param [Array] poslist [x,y]の配列群
  def initialize(poslist)
    init(poslist)
  end

  # 係数を計算して記憶
  # @param [Array] poslist [x,y]の配列群
  def init(poslist)
    @n = poslist.length - 1
    @h = []
    @b = []
    @d = []
    @g = []
    @u = []
    @q = []
    @s = []
    @r = []
    @x = []
    @y = []

    poslist.each do |px, py|
      @x.push(px.to_f)
      @y.push(py.to_f)
    end

    # step 1
    0.step(@n - 1, 1) do |i|
      @h[i] = @x[i + 1] - @x[i]
    end

    1.step(@n - 1, 1) do |i|
      @b[i] = 2.0 * (@h[i] + @h[i - 1])
      @d[i] = 3.0 * ((@y[i+1] - @y[i]) / @h[i] - (@y[i] - @y[i-1]) / @h[i-1])
    end

    # step 2
    @g[1] = @h[1] / @b[1]
    2.step(@n - 2, 1) do |i|
      @g[i] = @h[i] / (@b[i] - @h[i-1] * @g[i-1])
    end

    @u[1] = @d[1] / @b[1]
    2.step(@n - 1, 1) do |i|
      @u[i] = (@d[i] - @h[i-1] * @u[i-1]) / (@b[i] - @h[i-1] * @g[i-1])
    end

    # step 3
    @r[0] = 0.0
    @r[@n] = 0.0
    @r[@n-1] = @u[@n-1]
    (@n - 2).step(1, -1) do |i|
      @r[i] = @u[i] - @g[i] * @r[i+1]
    end

    # step 4
    0.step(@n - 1, 1) do |i|
      @q[i] = (@y[i+1] - @y[i]) / @h[i] - @h[i] * (@r[i+1] + 2.0 * @r[i]) / 3.0
      @s[i] = (@r[i+1] - @r[i]) / (3.0 * @h[i])
    end
  end

  # 補間値の計算
  # @param [Number] x1 補間値を求める値
  # @return [Number] 補間値
  def get_value(x1)
    i = -1.0

    i1 = 1
    while i1 < @n and i < 0.0
      i = i1 - 1 if x1 < @x[i1]
      i1 += 1
    end

    i = @n - 1 if i < 0.0

    xx = x1 - @x[i]
    y1 = @y[i] + xx * (@q[i] + xx * (@r[i] + @s[i] * xx))

    return y1
  end
end


if $0 == __FILE__
  # ----------------------------------------
  # 動作テスト

  # 四角を描画
  # @param [Number] x x座標
  # @param [Number] y y座標
  # @param [Number] sz サイズ
  # @param [Array] col 色。[r, g, b]の配列
  def draw_dot(x, y, sz, col)
    sz /= 2.0
    Window.drawBoxFill(x - sz, y - sz, x + sz, y + sz, col)
  end

  # 点を描画
  # @param [Number] x x座標
  # @param [Number] y y座標
  # @param [Array] col 色。[r, g, b]の配列
  def draw_dot2(x, y, color)
    Window.drawLine(x, y, x, y, color)
  end

  # 座標群を乱数で生成
  # @return [Array] [x0,y0],[x1,y1]...[xn,yn]の配列を生成して返す
  def get_poslist
    poslist = []
    w = Window.width
    h = Window.height
    0.step(w, 50) do |x|
      y = rand(h / 2) + (h / 4)
      poslist.push([x, y])
    end
    return poslist
  end

  fnt = Font.new(12)

  poslist = get_poslist
  spline = CubicSpline.new(poslist)

  # メインループ
  Window.loop do
    break if Input.keyPush?(K_ESCAPE) # ESCキーで終了

    if Input.mousePush?(M_LBUTTON)
      # マウスボタンが押されたら座標群を生成し直し
      poslist = get_poslist
      spline.init(poslist)
    end

    # 制御点を描画
    poslist.each do |x, y|
      draw_dot(x, y, 6, C_RED)
    end

    # 3次スプラインを描画
    end_x1 = poslist.last[0]
    x1 = 0.0
    x1_add = end_x1.to_f / (poslist.length * 16)
    while x1 <= (end_x1 + 0.01)
      y1 = spline.get_value(x1)
      draw_dot2(x1, y1, C_WHITE)
      x1 += x1_add
    end

    Window.drawFont(4, 4, "Please Click", fnt)
  end

end