スプライン補間をC言語で書いてみた

## はじめに
補間法の一つスプライン補間をC言語で書いてみました。
「何をいまさら」ですが、まあ見てください。
奥村晴彦先生の「C言語による標準アルゴリズム事典」を参考にしました。

## 方針
spline関数は表を読み、データーの繋目の条件をあらかじめ計算して置き、あとから呼び出されたときにその値を使って補間値を求めます。奥村先生の本では、一つのファイルに表を作る関数と、計算する関数を同居させました。
この本はアルゴリズムを説明するためにこうしてますが、IoT等組み込み分野では、この実装は不利です。

- 解析した結果は、初めのデーターとともにROMに配置したい
- 表の作成は、組み込み機ではなく、開発機で行いたい
- 複数の表を、交互に使いたい

これらを実現できるように書いてみます。
表（曲率データー）の作成と、実機での計算を分離したのが最大の特徴でが、同じマシン上でももちろん結構です。

## 表を作成するプログラム

``` makespline.c
/*
 * makespline.c:
 *      Make spline table
 *              Coded by fanfan (S.Makino)
 *              Copyright (c) 2021 Agri-hitech LLC.
 */

typedef struct {
        double in;
        double out;
        double diff;
} TBL;

// ここに表を書く
TBL data[] = {
        {0,0,0},
        {1,1,0},
        {2,2,0},
        {3,4,0},
        {4,8,0},
        {5,8,0},
        {6,4,0},
        {7,2,0},
        {8,1,0},
        {9,0,0},
};

/***********************************************************
        maketable.c -- スプライン補間データー作成
***********************************************************/

void maketable(TBL *T, int size)
{
        int i;
        double t;
        double h[size], d[size];

        T[0].diff = 0;  T[size - 1].diff = 0;  /* 両端点での y''(x) / 6 */
    for (i = 0; i < size - 1; i++) {
                h[i    ] =  T[i + 1].in - T[i].in;
                d[i + 1] = (T[i + 1].out - T[i].out) / h[i];
        }
        T[1].diff = d[2] - d[1] - h[0] * T[0].diff;
        d[1] = 2 * (T[2].in - T[0].in);
        for (i = 1; i < size - 2; i++) {
                t = h[i] / d[i];
                T[i + 1].diff = d[i + 2] - d[i + 1] - T[i].diff * t;
                d[i + 1] = 2 * (T[i + 2].in - T[i].in) - h[i] * t;
        }
        T[size - 2].diff -= h[size - 2] * T[size - 1].diff;
        for (i = size - 2; i > 0; i--)
                T[i].diff = (T[i].diff - h[i] * T[i + 1].diff) / d[i];
}
#define Nd (sizeof(data)/sizeof(TBL))

#include <stdio.h>
int main(){
        int i;

        maketable((void *)data,Nd);
        printf("// Auto generated by mkspline.c\n");
        printf("// DON'T Edit !!\n");
        printf("const TBL data[] = {\n");
        for (i = 0; i < Nd;i++) {
                printf("\t{%.16g,%.16g,%.16g},\n",
                                data[i].in,data[i].out,data[i].diff);
        }
        printf("};\n");

        return 0;
}
```
// ここに表を書く
TBL data[] = {
以下にに、データー点を書いてください。

{入力, 出力, 0},
と列挙します。
等間隔でなくても、入力が昇順になっていれば結構です。
３番目はこのプログラムが踏みつぶしますので、普通0を書いておきます。

## makesplineを実行する

``` sh
gcc makespline.c -o makespline
./makespline > data.c
```
data.cを見てみましょう。

``` data.c
// Auto generated by makespline.c
// DON'T Edit !!
const TBL data[] = {
        {0,0,0},
        {1,1,-0.0188679},
        {2,2,0.0754717},
        {3,4,0.716981},
        {4,8,-0.943396},
        {5,8,-0.943396},
        {6,4,0.716981},
        {7,2,0.0754717},
        {8,1,-0.0188679},
        {9,0,0},
};
```
プログラムに書いた表が変換されてます。
constがついていますので、ROMに配置され、3番目のデーターに計算結果が入ってます。

## スプライン本体

スプライン関数本体はこの表を取り込みます。
見てみましょう。

``` spline.c
/*
 * spline.c:
 *      Spline interporator
 *              Coded by fanfan (S.Makino)
 *              Copyright (c) 2021 Agri-hitech LLC.
 */
#define REAL float
//#define REAL double

typedef struct {
        REAL in;
        REAL out;
        REAL d;
} TBL;

#include "data.c"

/***********************************************************
        spline.c -- スプライン補間
***********************************************************/

#define  Nd sizeof(data)/sizeof(TBL)

float spline(float t, const TBL * T, int size)
{
        int i, j, k;
        REAL d, h;

        i = 0;  j = size - 1;
        while (i < j) {
                k = (i + j) / 2;
                if (T[k].in < t) i = k + 1;  else j = k;
        }
        if (i > 0) i--;
        h = T[i + 1].in - T[i].in;
        d = t - T[i].in;
        return (((T[i + 1].d - T[i].d) * d / h + T[i].d * 3) * d
                + ((T[i + 1].out - T[i].out) / h
                - (T[i].d * 2 + T[i + 1].d) * h)) * d + T[i].out;
}
#define DEBUG_MAIN
#ifdef DEBUG_MAIN
#include <stdio.h>
int main(){
        int i;
        REAL a;

        for (i = 0; i <= 90; i++) {
                a = (float) i / 10.0;
                printf("%.16g,%.16g\n",a, spline(a, data, Nd));
        }
}
#endif
```

`#define REAL float と宣言してますが、RAM/ROMに余裕があるなら、
`#define REAL doubleに変えてもいいです。精度が上がりますが、普通のセンサーシステムでは不要な範囲です。

実行してみましょう。

``` sh
gcc spline.c -o spline
./spline
```
結果をexcl等でグラフ化します。![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/245853/17bdab42-e17b-41c8-31d3-6141ebfbfedc.png)

## まとめ

- 目標は達成されました
- 構造体の使い方の例になっていると思います
- patchをくださった方ありがとうございます。

## 参考

https://qiita.com/YudaiSadakuni/items/e886ad237a1171e069a1