processing でスプラインを書いてみる

こんにちは。

前回、行列式の計算をするプログラムをprocessingで書いたので、これを利用してスプライン補間するプログラムを書いてみました。

今回、書いたのは、プロットした各点間を三次関数で補完するような曲線です。

こちらのサイトを参考にさせていただきました。 3 スプライン補間

コード

//processing

float[] X = {};
float[] Y = {};
int clickFlag = 1;

void setup()
{
  size(1000, 1000);
  fill(0, 0, 0);
}

void draw()
{
  if (mousePressed == false)
  {
    clickFlag = 1;
  }
  if (mousePressed == true && clickFlag == 1)
  {
    screenReset();
    X = append(X, mouseX);
    Y = append(Y, mouseY);
    for(int i = 0; i < X.length; i++)
    {
      ellipse(X[i], Y[i], 10, 10);
    }
    clickFlag = 0;
    spline(X, Y);
  }
}

//スクリーンを灰色にリセットする関数
void screenReset()
{
  fill(200, 200, 200);
  rect(0, 0, width, height);
}

//スプライン補間関数
void spline(float x[], float y[])
{
  int N = x.length;
  float[][] A = new float[N-1][N-1];
  float[][] B = new float[N-1][N-1];
  float[] h = new float[N];
  float[] v = new float[N];
  float[] u = new float[N];
  float[] a = new float[N];
  float[] b = new float[N];
  float[] c = new float[N];
  float[] d = new float[N];
  for (int i = 0; i < N-1; i++)
  {
    h[i] = x[i+1] - x[i];
  }
  for (int i = 1; i < N-1; i++)
  {
    v[i] = 6*((y[i+1] - y[i])/h[i] - (y[i]-y[i-1])/h[i-1]);
  }
  v[0] = 0;
  v[N-1] = 0;
  for (int i = 0; i < N-1; i++)
  {
    A[i][i] = 2*(h[i] + h[i+1]);
    if(i < N-2)
    {
      A[i][i+1] = h[i+1];
      A[i+1][i] = h[i+1];
    }
  }
  float X = determine(A);
  
  for (int ans = 0; ans < N-1; ans++)
  {
    for (int i = 0; i < N-1; i++)
    {
      for (int j = 0; j < N-1; j++)
      {
        B[i][j] = A[i][j];
      }
    }
    for(int i = 0; i < N-1; i++)
    {
      B[i][ans] = v[i+1];
    }
    u[ans+1] = determine(B)/X;
  } 
  u[0] = 0;
  u[N-1] = 0;

  for (int i = 0; i < N-1; i++)
  {
    a[i] = (u[i+1]-u[i])/(6*h[i]);
    b[i] = u[i]/2;
    c[i] = (y[i+1]-y[i])/h[i]-h[i]*(2*u[i]+u[i+1])/6;
    d[i] = y[i];
  }
  stroke(0, 0, 0);
  for(int i = 0; i < x[N-1]; i++)
  {
    for(int k = 1; k < N; k++)
    {
      if(x[k-1] < i && i < x[k])
      {
        line(i, (int)(a[k-1]*pow((i-x[k-1]), 3) + b[k-1]*pow((i-x[k-1]), 2) + c[k-1]*(i-x[k-1]) + d[k-1]), i+1, (int)(a[k-1]*pow((i+1-x[k-1]), 3) + b[k-1]*pow((i+1-x[k-1]), 2) + c[k-1]*(i+1-x[k-1]) + d[k-1]));
      }
    }
  }
}

//行列式の計算をする関数
float determine(float m[][])
{
    float result = 0;
    int n = m.length;
    if(n==1)
    {
      result = m[0][0];
    }
    else if(n==2)
    {
      result = m[0][0]*m[1][1]-m[0][1]*m[1][0];
    }
    else if(n==3)
    {
      result =  m[0][0]*m[1][1]*m[2][2]-m[0][0]*m[1][2]*m[2][1]
            +m[0][1]*m[1][2]*m[2][0]-m[0][1]*m[1][0]*m[2][2]
            +m[0][2]*m[1][0]*m[2][1]-m[0][2]*m[1][1]*m[2][0];
    }
    if(n > 3)
    {
      float[][][] _m = new float[n][n-1][n-1];
      float[] c = new float[n];
      for(int k = 0; k < n; k++)
      {
        c[k] = m[k][0];
        for(int j = 1; j < n; j++)
        {
            for(int i = 0; i < k; i++)
            {
              _m[k][i][j-1] = m[i][j];
            }
            for(int i = k+1; i < n; i++)
            {
              _m[k][i-1][j-1] = m[i][j];
            }
        }
        result += c[k] * pow(-1, k) * determine(_m[k]);
      }
    }
    return result;
}

実行したら画面を左から右へとクリックしていってください。その点がプロットされ、また、補間されるはずです。

こんな感じです。

点の数が10個以上になると計算に時間がかかっていることが目に見えて分かりますね...

そもそも、こういった計算はライブラリに任せるのが一番だとは思いますが、１から書いて理解が深まったので良しとしましょう。