I used this method to delete a gradient effect from uneven lighting of some of my book scans.

Here is a simple method for curve fitting various polynomials. This works for non linear equations as long as the actual coeficients are linear.

For example you can NOT fit the following type of function:

Note that the coeficients is an argument to a nonlinear function. related link

The specific one done here is a bicubic polynomial expressed by the following.

We are going to fit the following curve:

Or more compactly:

import cern.colt.matrix.linalg.*; import cern.colt.matrix.*; import cern.colt.matrix.impl.*; import edu.umbc.cs.maple.utils.*; public DoubleMatrix2D calcFit(ArrayList<Point> measures) { int cols = 16; // number of coeficients int rows = measures.size(); // size of your dataset DenseDoubleMatrix2D M = new DenseDoubleMatrix2D(rows,cols); DenseDoubleMatrix2D P = new DenseDoubleMatrix2D(rows,1); for (int k=0;k<measures.size();k++) { Point measure = measures.get(k); // measures from your dataset P.set(k,0,measure.getValue()); // Point contains values P , x and y double x = measure.getX(); double y = measure.getY(); // populate the matrix with the appropriate polynomial values // from the bicubic above int col=0; for (int i=0;i<4;i++) for (int j=0;j<4;j++) { double cal = pow(x,i)*pow(y,j); M.set(k,col++,cal); } } // solve SVD SingularValueDecomposition svd = new SingularValueDecomposition(M); DoubleMatrix2D U = svd.getU(); DoubleMatrix2D V = svd.getV(); DoubleMatrix2D S = svd.getS(); // S is the diagonal matrix // Invert by replacing its diagonal elements by their reciprocals for (int i=0;i<16;i++) { double val = 1.0/S.get(i,i); if (S.get(i,i) == 0) // avoid numeric underflow val = 0.0; S.set(i,i, val); } Algebra alg = new Algebra(); // X pseudoinverse DoubleMatrix2D VS = alg.mult(V,S); DoubleMatrix2D Minv = alg.mult(VS,alg.transpose(U)); // matrix has coeficients DoubleMatrix2D C = alg.mult(Minv, P); // If you would like the closeness of the fit use the following DoubleMatrix2D eP = alg.mult(M,C); // estimated values ColtUtils utils = new ColtUtils(); DoubleMatrix2D dP = utils.minus(P,eP); // difference between estimated and observed values double sum = utils.dotproduct(utils.getcol(dP,0),utils.getcol(dP,0)); System.out.println("sum = " + sum); // how close was our fit? return(C); }

In practice you may need to normalize your measurements if they are large to avoid numeric under and/or over flow.