By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
 428,759 Members | 1,727 Online + Ask a Question
Need help? Post your question and get tips & solutions from a community of 428,759 IT Pros & Developers. It's quick & easy.

# best way(s) to fit a polynomial to points?

 P: n/a I'm trying to do some georeferencing - using points of known location (ground control points, GCPs) on an image to develop a polynomial that can be used to approximate the locations of other points. I usually use the Python bindings to GDAL http://www.remotesensing.org/gdal/ to manipulate geospatial data but there are no Python bindings for GDALCreateGCPTransformer(). http://remotesensing.org/cgi-bin/cvs...ype=text/plain I tried using it from C (both through GDAL and directly to the GRASS code) and failed. I'd rather have a Python-based solution anyway, so I'm searching for appropriate tools. The ones I'm finding (like matplotlib.mlab.polyfit() and PyBLD's polfit()) only deal with single variable equations (f(x) = X). I need something that operates on two variables (f(x, y) = X). (I'd have two polynomials; one would calculate the northing and the other would do the easting). I need to be able to fit at least third degree polynomials. The solution does not need to be especially fast. A pure Python solution is preferred but I'll take about anything that works right now. (I'm trying to show someone how to do this. I've only used proprietary software to do it but I've been meaning to write my own.) Any suggestions? Thank you. --kyler Jul 18 '05 #1
Share this Question
4 Replies

 P: n/a In article , Kyler Laird wrote:I'm trying to do some georeferencing - using points of knownlocation (ground control points, GCPs) on an image to developa polynomial that can be used to approximate the locations ofother points.I usually use the Python bindings to GDAL http://www.remotesensing.org/gdal/to manipulate geospatial data but there are no Python bindingsfor GDALCreateGCPTransformer(). http://remotesensing.org/cgi-bin/cvs...ype=text/plainI tried using it from C (both through GDAL and directly tothe GRASS code) and failed.I'd rather have a Python-based solution anyway, so I'msearching for appropriate tools. The ones I'm finding (likematplotlib.mlab.polyfit() and PyBLD's polfit()) only deal withsingle variable equations (f(x) = X). I need something thatoperates on two variables (f(x, y) = X). (I'd have twopolynomials; one would calculate the northing and the otherwould do the easting). I need to be able to fit at leastthird degree polynomials.The solution does not need to be especially fast. A purePython solution is preferred but I'll take about anythingthat works right now. (I'm trying to show someone how to dothis. I've only used proprietary software to do it but I'vebeen meaning to write my own.) Jul 18 '05 #2

 P: n/a cl****@lairds.com (Cameron Laird) writes: "I tried using it ... and failed." I naturally wonder *how*you (it, more likely) failed--the specific details of whatyou actually observed. I was generating the test points with a simple algorithm. I finally discovered that the fitting algorithm fails if the points are (close to being) colinear. I can imagine, though, that if youhad the time to understand and describe what happened withprecision, it'd be easier fixing it than talking about it. I'll append the code that I successfully used. Will the new (x,y)-s be interior to the convex hull formed bythe tabulated ones, or might they be outside? It's good practice to make them interior (have ground control points around the perimeter of the region of interest) but it isn't always the case. The solution should be general enough to work for outside points but it's understood that the error for them is likely to be high. Howsoon does this thing need to go live? It's something I meant to do during my Remote Sensing class last semester and then someone asked about it in another news group recently. It's not urgent. I just wanted to solve it. (I have some aerial photos that I georeferenced using ERDAS Imagine but I want to use my own code so that I can do more with them.) I'm happy with the C code, but I'd still prefer a Python solution. --kyler ================================================== ============ #include struct Control_Points { int count; double *e1; double *n1; double *e2; double *n2; int *status; }; /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS. THESE FUNCTIONS EXPECT SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */ struct MATRIX { int n; /* SIZE OF THIS MATRIX (N x N) */ double *v; }; /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */ #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1] #define MSUCCESS 1 /* SUCCESS */ #define MNPTERR 0 /* NOT ENOUGH POINTS */ #define MUNSOLVABLE -1 /* NOT SOLVABLE */ #define MMEMERR -2 /* NOT ENOUGH MEMORY */ #define MPARMERR -3 /* PARAMETER ERROR */ #define MINTERR -4 /* INTERNAL ERROR */ #define MAXORDER 3 #define CPLFree free #define CPLCalloc calloc /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */ #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1] #define MSUCCESS 1 /* SUCCESS */ #define MNPTERR 0 /* NOT ENOUGH POINTS */ #define MUNSOLVABLE -1 /* NOT SOLVABLE */ #define MMEMERR -2 /* NOT ENOUGH MEMORY */ #define MPARMERR -3 /* PARAMETER ERROR */ #define MINTERR -4 /* INTERNAL ERROR */ /************************************************** *************************/ /* FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS */ /************************************************** *************************/ static int calccoef(struct Control_Points *,double *,double *,int); static int calcls(struct Control_Points *,struct MATRIX *, double *,double *,double *,double *); static int exactdet(struct Control_Points *,struct MATRIX *, double *,double *,double *,double *); static int solvemat(struct MATRIX *,double *,double *,double *,double *); static double term(int,double,double); /************************************************** *************************/ /* TRANSFORM A SINGLE COORDINATE PAIR. */ /************************************************** *************************/ static int CRS_georef ( double e1, /* EASTINGS TO BE TRANSFORMED */ double n1, /* NORTHINGS TO BE TRANSFORMED */ double *e, /* EASTINGS TO BE TRANSFORMED */ double *n, /* NORTHINGS TO BE TRANSFORMED */ double E[], /* EASTING COEFFICIENTS */ double N[], /* NORTHING COEFFICIENTS */ int order /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE ORDER USED TO CALCULATE THE COEFFICIENTS */ ) { double e3, e2n, en2, n3, e2, en, n2; switch(order) { case 1: *e = E + E * e1 + E * n1; *n = N + N * e1 + N * n1; break; case 2: e2 = e1 * e1; n2 = n1 * n1; en = e1 * n1; *e = E + E * e1 + E * n1 + E * e2 + E * en + E * n2; *n = N + N * e1 + N * n1 + N * e2 + N * en + N * n2; break; case 3: e2 = e1 * e1; en = e1 * n1; n2 = n1 * n1; e3 = e1 * e2; e2n = e2 * n1; en2 = e1 * n2; n3 = n1 * n2; *e = E + E * e1 + E * n1 + E * e2 + E * en + E * n2 + E * e3 + E * e2n + E * en2 + E * n3; *n = N + N * e1 + N * n1 + N * e2 + N * en + N * n2 + N * e3 + N * e2n + N * en2 + N * n3; break; default: return(MPARMERR); break; } return(MSUCCESS); } /************************************************** *************************/ /* COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS */ /************************************************** *************************/ static int CRS_compute_georef_equations (struct Control_Points *cp, double E12[], double N12[], double E21[], double N21[], int order) { double *tempptr; int status; if(order < 1 || order > MAXORDER) return(MPARMERR); /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */ status = calccoef(cp,E12,N12,order); if(status != MSUCCESS) return(status); /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */ tempptr = cp->e1; cp->e1 = cp->e2; cp->e2 = tempptr; tempptr = cp->n1; cp->n1 = cp->n2; cp->n2 = tempptr; /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */ status = calccoef(cp,E21,N21,order); /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */ tempptr = cp->e1; cp->e1 = cp->e2; cp->e2 = tempptr; tempptr = cp->n1; cp->n1 = cp->n2; cp->n2 = tempptr; return(status); } /************************************************** *************************/ /* COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS */ /************************************************** *************************/ static int calccoef (struct Control_Points *cp, double E[], double N[], int order) { struct MATRIX m; double *a; double *b; int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */ int status, i; /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */ for(i = numactive = 0 ; i < cp->count ; i++) { if(cp->status[i] > 0) numactive++; } /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE A TRANSFORMATION OF THIS ORDER */ m.n = ((order + 1) * (order + 2)) / 2; if(numactive < m.n) return(MNPTERR); /* INITIALIZE MATRIX */ m.v = (double *)CPLCalloc(m.n*m.n,sizeof(double)); if(m.v == NULL) { return(MMEMERR); } a = (double *)CPLCalloc(m.n,sizeof(double)); if(a == NULL) { CPLFree((char *)m.v); return(MMEMERR); } b = (double *)CPLCalloc(m.n,sizeof(double)); if(b == NULL) { CPLFree((char *)m.v); CPLFree((char *)a); return(MMEMERR); } if(numactive == m.n) status = exactdet(cp,&m,a,b,E,N); else status = calcls(cp,&m,a,b,E,N); /* CPLFree((char *)m.v); CPLFree((char *)a); CPLFree((char *)b); */ return(status); } /************************************************** *************************/ /* CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. */ /************************************************** *************************/ static int exactdet ( struct Control_Points *cp, struct MATRIX *m, double a[], double b[], double E[], /* EASTING COEFFICIENTS */ double N[] /* NORTHING COEFFICIENTS */ ) { int pntnow, currow, j; currow = 1; for(pntnow = 0 ; pntnow < cp->count ; pntnow++) { if(cp->status[pntnow] > 0) { /* POPULATE MATRIX M */ for(j = 1 ; j <= m->n ; j++) { M(currow,j) = term(j,cp->e1[pntnow],cp->n1[pntnow]); } /* POPULATE MATRIX A AND B */ a[currow-1] = cp->e2[pntnow]; b[currow-1] = cp->n2[pntnow]; currow++; } } if(currow - 1 != m->n) { return(MINTERR); } return(solvemat(m,a,b,E,N)); } /************************************************** *************************/ /* CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. THIS ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS. */ /************************************************** *************************/ static int calcls ( struct Control_Points *cp, struct MATRIX *m, double a[], double b[], double E[], /* EASTING COEFFICIENTS */ double N[] /* NORTHING COEFFICIENTS */ ) { int i, j, n, numactive = 0; /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */ for(i = 1 ; i <= m->n ; i++) { for(j = i ; j <= m->n ; j++) M(i,j) = 0.0; a[i-1] = b[i-1] = 0.0; } /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */ for(n = 0 ; n < cp->count ; n++) { if(cp->status[n] > 0) { numactive++; for(i = 1 ; i <= m->n ; i++) { for(j = i ; j <= m->n ; j++) M(i,j) += term(i,cp->e1[n],cp->n1[n]) * term(j,cp->e1[n],cp->n1[n]); a[i-1] += cp->e2[n] * term(i,cp->e1[n],cp->n1[n]); b[i-1] += cp->n2[n] * term(i,cp->e1[n],cp->n1[n]); } } } if(numactive <= m->n) return(MINTERR); /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */ for(i = 2 ; i <= m->n ; i++) { for(j = 1 ; j < i ; j++) M(i,j) = M(j,i); } return(solvemat(m,a,b,E,N)); } /************************************************** *************************/ /* CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER ORDER\TERM 1 2 3 4 5 6 7 8 9 10 1 e0n0 e1n0 e0n1 2 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 3 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3 */ /************************************************** *************************/ static double term (int term, double e, double n) { switch(term) { case 1: return((double)1.0); case 2: return((double)e); case 3: return((double)n); case 4: return((double)(e*e)); case 5: return((double)(e*n)); case 6: return((double)(n*n)); case 7: return((double)(e*e*e)); case 8: return((double)(e*e*n)); case 9: return((double)(e*n*n)); case 10: return((double)(n*n*n)); } return((double)0.0); } /************************************************** *************************/ /* SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED GAUSSIAN ELIMINATION METHOD. | M11 M12 ... M1n | | E0 | | a0 | | M21 M22 ... M2n | | E1 | = | a1 | | . . . . | | . | | . | | Mn1 Mn2 ... Mnn | | En-1 | | an-1 | and | M11 M12 ... M1n | | N0 | | b0 | | M21 M22 ... M2n | | N1 | = | b1 | | . . . . | | . | | . | | Mn1 Mn2 ... Mnn | | Nn-1 | | bn-1 | */ /************************************************** *************************/ static int solvemat (struct MATRIX *m, double a[], double b[], double E[], double N[]) { int i, j, i2, j2, imark; double factor, temp; double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */ for(i = 1 ; i <= m->n ; i++) { j = i; /* find row with largest magnitude value for pivot value */ pivot = M(i,j); imark = i; for(i2 = i + 1 ; i2 <= m->n ; i2++) { temp = fabs(M(i2,j)); if(temp > fabs(pivot)) { pivot = M(i2,j); imark = i2; } } /* if the pivot is very small then the points are nearly co-linear */ /* co-linear points result in an undefined matrix, and nearly */ /* co-linear points results in a solution with rounding error */ if(pivot == 0.0) { return(MUNSOLVABLE); } /* if row with highest pivot is not the current row, switch them */ if(imark != i) { for(j2 = 1 ; j2 <= m->n ; j2++) { temp = M(imark,j2); M(imark,j2) = M(i,j2); M(i,j2) = temp; } temp = a[imark-1]; a[imark-1] = a[i-1]; a[i-1] = temp; temp = b[imark-1]; b[imark-1] = b[i-1]; b[i-1] = temp; } /* compute zeros above and below the pivot, and compute values for the rest of the row as well */ for(i2 = 1 ; i2 <= m->n ; i2++) { if(i2 != i) { factor = M(i2,j) / pivot; for(j2 = j ; j2 <= m->n ; j2++) M(i2,j2) -= factor * M(i,j2); a[i2-1] -= factor * a[i-1]; b[i2-1] -= factor * b[i-1]; } } } /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */ for(i = 1 ; i <= m->n ; i++) { E[i-1] = a[i-1] / M(i,i); N[i-1] = b[i-1] / M(i,i); } return(MSUCCESS); } main() { int i; double E12, N12, E21, N21; int order = 1; int ret; double x, y, X, Y; struct Control_Points CPs; CPs.count = 10; CPs.e1 = (double *)malloc(sizeof(double) * CPs.count); CPs.n1 = (double *)malloc(sizeof(double) * CPs.count); CPs.e2 = (double *)malloc(sizeof(double) * CPs.count); CPs.n2 = (double *)malloc(sizeof(double) * CPs.count); CPs.status = (int *)malloc(sizeof(int) * CPs.count); i = 0; /* Set some GCPs. */ CPs.e1[i] = 0; CPs.n1[i] = 0; CPs.e2[i] = 0; CPs.n2[i] = 0; CPs.status[i] = 1; i++; CPs.e1[i] = 10; CPs.n1[i] = 30; CPs.e2[i] = 100; CPs.n2[i] = 200; CPs.status[i] = 1; i++; CPs.e1[i] = 15; CPs.n1[i] = 35; CPs.e2[i] = 150; CPs.n2[i] = 220; CPs.status[i] = 1; i++; for(i = 0; i < 30; i++) { E12[i] = N12[i] = E21[i] = N21[i] = 0; } CRS_compute_georef_equations( &CPs, E12, N12, E21, N21, order ); for(i = 0; i < CPs.count; i++) { x = i * 100; y = i * 200; CRS_georef(x, y, &X, &Y, E21, N21, order); printf("(%lf, %lf) -> (%lf, %lf)\n", x, y, X, Y); CRS_georef(X, Y, &x, &y, E12, N12, order); printf("(%lf, %lf) -> (%lf, %lf)\n", X, Y, x, y); } } Jul 18 '05 #3

 P: n/a In article <3o************@jowls.lairds.org>, Kyler Laird wrote: Jul 18 '05 #4

 P: n/a Hi, How about using `scipy.optimize` at http://www.scipy.org . From scipy_tutorial.pdf, it saysA collection of general-purpose optimization routines. fmin -- Nelder-Mead Simplex algorithm (uses only function calls) fmin_powell -- Powells (modified) level set method (uses onlyfunction calls) fmin_bfgs -- Quasi-Newton method (can use function and gradient) fmin_ncg -- Line-search Newton Conjugate Gradient (can usefunction, gradient and hessian). leastsq -- Minimize the sum of squares of M equations in N unknowns given a starting estimate. Note that if you want Win binary, you'll find them under http://www.scipy.org/site_content/downloads where you cannot reach by clicking links, I suppose. ----- Natsu Jul 18 '05 #5

### This discussion thread is closed

Replies have been disabled for this discussion. 