スプライン補間

以下の１０個のデータの点を曲線でつなぎグラフを作るとします。

 1,1,1,1,2,2,2,3,3,7

さて、3 と 7 の真ん中の数字はいくつになるでしょう？

直線で結ぶのであれば 3 と 7 を足して 2 で割った 5 になります。
しかし、曲線で結ぶわけですから直線でやる場合とは同じようにいきません。
この答えは 3 より以前の 1 や 2 といった数字の並びが影響してきます。
おそらく 5 よりは小さい値であることは想像できます。
この間の数値を計算で求めるのがスプライン補間です。
スプライン補間を使うとこの値がだいたい 4.68 であることがわかります。

スプライン補間はどういう場合に役立つでしょうか？
例えば、データ数１３個のＷＡＶＥデータがあるとします。
このデータ数を１３から１６に増やしたいと思うとき。
１３個のデータの端から端まで１６分割してその間の数値を求めれば、
ＷＡＶＥの波形を崩さずに１６個のデータに変換できます。
別の例では、例えば横幅１００ピクセルの画像データがあったとします。
これを横幅１１０ピクセルの画像データに引き伸ばしたいとします。
この場合も１００ピクセルのデータの端から端を１１０分割して、
間のデータを求めれば１１０ピクセルでの各ピクセルでのデータ値がわかります。
こういった”間の数値”の予想は計算で行うことが可能です。
引き伸ばしだけではなく縮小についてもスプライン補間が利用できそうです。

スプライン補間のアルゴリズムは各データの x 値の間隔が任意でも計算できますが、
下のサンプルは x 値の間隔を 1.0 で固定したスプライン補間のコードです。
10 個のデータを使用して 100 個の滑らかなデータに引き伸ばしています。
結果を Excel などでグラフにしてみると効果がわかりやすいと思います。

/*
*  スプライン補間のサンプル
*  (ただし x 軸の間隔は 1.0 で固定である。)
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
*  変換前の配列 b と配列数 num_b
*  変換後の配列 a と配列数 num_a
*  b と a は関数呼び出し前に実体を用意すること
*/
int spline_interpolation(float *b, int num_b, float *a, int num_a){
  int i;
  int num1 = num_b - 1;
  int num2 = num_a - 1;

  float *d = NULL;
  float *g = NULL;
  float *u = NULL;
  float *r = NULL;

  d = (float *)malloc(num1 * sizeof(float));
  g = (float *)malloc(num1 * sizeof(float));
  u = (float *)malloc(num1 * sizeof(float));
  r = (float *)malloc((num1 + 1) * sizeof(float));

  for(i = 1; i < num1; i++){
    d[i] = 3.0f * (b[i + 1] + b[i - 1] - 2.0f * b[i]);
  }

  g[1] = 1.0f / 4.0f;
  
  for(i = 2; i < (num1 - 1); i++){
    g[i] = 1.0f / (4.0f - g[i - 1]); 
  }

  u[1] = d[1] / 4.0f;
  
  for(i = 2; i < num1; i++){
    u[i] = (d[i] - u[i - 1]) / (4.0f - g[i - 1]);
  }

  r[num1] = 0.0f;
  r[num1 - 1] = u[num1 - 1];
  
  for(i = num1 - 2; i > 0; i--){
    r[i] = u[i] - g[i] * r[i + 1]; 
  }

  r[0] = 0.0f;

  for(i = 0; i < num2; i++){
    double dd = (double)num1 / (double)num2 * (double)i;
    double sei;
    double syo = modf(dd, &sei);
    int t = (int)sei;
    float v = (float)syo;

    float q = (b[t + 1] - b[t]) - (r[t + 1] + 2.0f * r[t]) / 3.0f;
    float s = (r[t + 1] - r[t]) / 3.0f;

    a[i] = b[t] + v * (q + v * (r[t] + s * v));
  }

  a[num2] = b[num1];

  free(d);
  free(g);
  free(u);
  free(r);

  return 1;
}

/*
*   要素 10 の配列をスプライン補間にて要素 100 の配列に展開
*/
int main(){
  float a[10] = {0.0f, 1.0f, 4.0f, 9.0f, 4.0f, 3.0f, 4.0f, 4.0f, 3.0f, 1.0f};
  float b[100];
  int i;

  spline_interpolation(a, 10, b, 100);

  for(i = 0 ; i < 100; i++){
    printf("%03d %f\n", i, b[i]);
  }

  return 0;
}