From 38ba969de8327d100130a5e77bce3e6cd4cc5d63 Mon Sep 17 00:00:00 2001 From: Thierry Leconte Date: Tue, 11 Nov 2003 20:56:01 +0000 Subject: [PATCH] use polynomial regression insteed of linear for calibration --- image.c | 32 ++++------- reg.c | 168 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 180 insertions(+), 20 deletions(-) create mode 100644 reg.c diff --git a/image.c b/image.c index 5aeb713..bc5ed36 100644 --- a/image.c +++ b/image.c @@ -25,38 +25,30 @@ #include #include +#define REGORDER 3 typedef struct { -double slope; -double offset; +double cf[REGORDER+1] ; } rgparam; static void rgcomp(double x[16], rgparam *rgpr) { /* 0.106,0.215,0.324,0.434,0.542,0.652,0.78,0.87 ,0.0 */ const double y[9] = { 31.1,63.0,95.0,127.2,158.9,191.1,228.6,255.0, 0.0 }; -const double yavg=(y[0]+y[1]+y[2]+y[4]+y[5]+y[6]+y[7]+y[8])/9.0; -double xavg; -double sxx,sxy; -int i; - -for(i=0,xavg=0.0;i<9;i++) - xavg+=x[i]; -xavg/=9; -for(i=0,sxx=0.0;i<9;i++) { - float t=x[i]-xavg; - sxx+=t*t; -} -for(i=0,sxy=0.0;i<9;i++) { - sxy+=(x[i]-xavg)*(y[i]-yavg); -} -rgpr->slope=sxy/sxx; -rgpr->offset=yavg-rgpr->slope*xavg; +extern void polyreg(int m,int n,double x[],double y[],double c[]); +polyreg(REGORDER,9,x,y,rgpr->cf); } static double rgcal(float x,rgparam *rgpr) { -return(rgpr->slope*x+rgpr->offset); +double y,p; +int i; + +for(i=0,y=0.0,p=1.0;icf[i]*p; + p=p*x; +} +return(y); } diff --git a/reg.c b/reg.c new file mode 100644 index 0000000..40a2368 --- /dev/null +++ b/reg.c @@ -0,0 +1,168 @@ +/* --------------------------------------------------------------------------- + + Polynomial regression, freely adapted from : + + NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995 + Algorithm translated to C by: Dr. Norman Fahrer + NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 + Prentice Hall, International Editions: ISBN 0-13-625047-5 + This free software is compliments of the author. + E-mail address: in%"mathews@fullerton.edu" +*/ + +#include + +#define DMAX 5 /* Maximum degree of polynomial */ +#define NMAX 10 /* Maximum number of points */ + +static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]); + +void polyreg(int M, int N, double X[], double Y[], double C[]) +{ + int R, K, J; /* Loop counters */ + double A[DMAX][DMAX]; /* A */ + double B[DMAX]; + double P[2*DMAX+1]; + double x, y; + double p; + + /* Zero the array */ + + for (R = 0; R < M+1; R++) B[R] = 0; + + /* Compute the column vector */ + + for (K = 0; K < N; K++) + { + y = Y[K]; + x = X[K]; + p = 1.0; + + for( R = 0; R < M+1; R++ ) + { + B[R] += y * p; + p = p*x; + } + } + + /* Zero the array */ + + for (J = 1; J <= 2*M; J++) P[J] = 0; + + P[0] = N; + + /* Compute the sum of powers of x_(K-1) */ + + for (K = 0; K < N; K++) + { + x = X[K]; + p = X[K]; + + for (J = 1; J <= 2*M; J++) + { + P[J] += p; + p = p * x; + } + } + + /* Determine the matrix entries */ + + for (R = 0; R < M+1; R++) + { + for( K = 0; K < M+1; K++) A[R][K] = P[R+K]; + } + + /* Solve the linear system of M + 1 equations : A*C = B + for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */ + FactPiv(M+1, A, B, C); + +} /* end main */ + +/*--------------------------------------------------------*/ + +static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) +{ + + int K, P, C, J; /* Loop counters */ + int Row[NMAX]; /* Field with row-number */ + double X[DMAX], Y[DMAX]; + double SUM, DET = 1.0; + int T; + + /* Initialize the pointer vector */ + + for (J = 0; J< N; J++) Row[J] = J; + + /* Start LU factorization */ + + for (P = 0; P < N - 1; P++) { + + /* Find pivot element */ + + for (K = P + 1; K < N; K++) { + if ( fabs(A[Row[K]][P]) > fabs(A[Row[P]][P]) ) { + /* Switch the index for the p-1 th pivot row if necessary */ + T = Row[P]; + Row[P] = Row[K]; + Row[K] = T; + DET = - DET; + } + + } /* End of simulated row interchange */ + + if (A[Row[P]][P] == 0) { + printf("The matrix is SINGULAR !\n"); + printf("Cannot use algorithm --> exit\n"); + exit(1); + } + + /* Multiply the diagonal elements */ + + DET = DET * A[Row[P]][P]; + + /* Form multiplier */ + + for (K = P + 1; K < N; K++) { + A[Row[K]][P]= A[Row[K]][P] / A[Row[P]][P]; + + /* Eliminate X_(p-1) */ + + for (C = P + 1; C < N + 1; C++) { + A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C]; + } + } + + } /* End of L*U factorization routine */ + + DET = DET * A[Row[N-1]][N-1]; + + /* Start the forward substitution */ + + for(K = 0; K < N; K++) Y[K] = B[K]; + + Y[0] = B[Row[0]]; + for ( K = 1; K < N; K++) { + SUM =0; + for ( C = 0; C <= K -1; C++) SUM += A[Row[K]][C] * Y[C]; + Y[K] = B[Row[K]] - SUM; + } + + /* Start the back substitution */ + + X[N-1] = Y[N-1] / A[Row[N-1]][N-1]; + + for (K = N - 2; K >= 0; K--) { + SUM = 0; + for (C = K + 1; C < N; C++) { + SUM += A[Row[K]][C] * X[C]; + } + + X[K] = ( Y[K] - SUM) / A[Row[K]][K]; + + } /* End of back substitution */ + + /* Output */ + for( K = 0; K < N; K++) + Cf[K]=X[K]; +} +