距離校正はつらいよ（6）

こんにちは、ソフト班Shinpukuです。

ここを参考にして作った「n次多項式近似のプログラム」です。

const int num = 16;

double distance[num] = {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5};

double pixels[num] = { 21, 27, 31, 37, 42, 44, 46, 50, 62, 68, 74, 76, 80, 82, 84, 84, 86};

// n次多項式近似 ///////////////////////////////////////////////////////
const int N = 5; // N-1次の多項式近似
double x[N];
double a[N][N+1];
double b[1][N+1];

for (int i = 0; i < N; i++) {
    for (int j = 0; j < N + 1; j++) {
        a[i][j] = 0.0;
    }
}
for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
        for (int k = 0; k < num; k++) {
            a[i][j] += pow((double)pixels[k], (double)(i + j));
        }
    }
}
for (int i = 0; i < N; i++) {
    for (int j = 0; j < num; j++) {
        a[i][N] += pow((double)pixels[j], (double)i) * distance[j];
    }
}
for (int i = 0; i < N; i++) {
    double m = 0;
    int pivot = i;
    for (int j = i; j < N; j++) {
        if (fabs(a[j][i]) > m) {
            m = fabs(a[j][i]);
            pivot = j;
        }
    }
    if (pivot != i) {
        for(int j = 0; j < N+1; j++) {
            b[0][j] = a[i][j];
            a[i][j] = a[pivot][j];
            a[pivot][j] = b[0][j];
        }
    }
}
for (int k = 0; k < N; k++) {
    double p = a[k][k];
    a[k][k] = 1;
    for (int j = k+1; j < N+1; j++) a[k][j] /= p;
    for (int i = k+1; i < N; i++) {
        double q = a[i][k];
        for (int j = k+1; j < N+1; j++) a[i][j] -= q*a[k][j];
        a[i][k] = 0;
    }
}
for (int i = N-1; i >= 0; i--) {
    x[i] = a[i][N];
    for (int j = N-1; j > i; j--) x[i] -= a[i][j]*x[j];
}
//////////////////////////////////////////////////////////////////////

これで、
r = 中心からのピクセル数
y = x[4]*r*r*r*r+ x[3]*r*r*r+ x[2]*r*r+ x[1]*r+ x[0];
で距離yがでます。