linear regression example w/ GSL

Ok so this is kind of random, but I’ve been playing with some functionality that I thought people might like to know about, which is doing linear regression to get the best parameters for a mapping equation. In my case, it’s for mapping tracked coordinates “eye” to “screen” coordinates. The mapping is non linear, and I tried warping which didn’t work so well. this uses GSL (http://www.gnu.org/software/gsl/) gnu scientific library but there might be other simpler libraries (this one worked quite fast). The linear regression allows you to find the best fit for any equation, ie:

Ax + Bx^2 + Cy + Dy^2 + E

it finds the A,B,C,D,E that best fit the mapping, and you need to do it for X and Y.

I’ll post a mac version that works (with GSL) shortly, but I thought I’d put this up cause it’s pretty cool:

in .h:

  
  
void calculateWeights(vector <ofPoint> eyePoints, vector <ofPoint>  screenPoints);  
	  
vector <ofPoint> follow;  
vector <ofPoint> myMouse;  
	  
double	cxfit[8];  
double	cyfit[8];  
bool	bBeenFit;  
  

in .cpp

relevant includes:

  
  
#include <gsl_math.h>  
#include <gsl_test.h>  
#include <gsl_fit.h>  
#include <gsl_multifit.h>  
#include <gsl_multifit_nlin.h>  
#include <gsl_blas.h>  
#include <gsl_ieee_utils.h>  
  

fitting function:

  
  
void testApp::calculateWeights(vector <ofPoint> eyePoints, vector <ofPoint>  screenPoints){  
	  
	int length = eyePoints.size();  
	  
	int nTerms = 8;  
	  
	gsl_matrix * x = gsl_matrix_alloc(length,nTerms);  
	gsl_vector * yx = gsl_vector_alloc(length);  
	gsl_vector * yy = gsl_vector_alloc(length);  
	gsl_vector * w = gsl_vector_alloc(nTerms);  
	  
	double * ptr;  
	double * ptrScreenX;  
	double * ptrScreenY;  
  
	  
	ptr = gsl_matrix_ptr(x,0,0);  
	ptrScreenX = gsl_vector_ptr(yx,0);  
	ptrScreenY = gsl_vector_ptr(yy,0);  
	  
	  
	for (int i = 0; i < length; i++){  
  
		float xPosEye = eyePoints[i].x;  
		float yPosEye = eyePoints[i].y;  
		*ptr++ = xPosEye;  
		*ptr++ = xPosEye*xPosEye;  
		*ptr++ = yPosEye;  
		*ptr++ = yPosEye*yPosEye;  
		*ptr++ = xPosEye*yPosEye;  
		*ptr++ = xPosEye*xPosEye*xPosEye;  
		*ptr++ = yPosEye*yPosEye*yPosEye;  
		*ptr++ = 1;  
		  
		*ptrScreenX++ = screenPoints[i].x;  
		*ptrScreenY++ = screenPoints[i].y;  
		  
	}  
	  
	  
	gsl_vector *cx = gsl_vector_calloc(nTerms);  
	gsl_vector *cy = gsl_vector_calloc(nTerms);  
	  
	  
    gsl_matrix *cov = gsl_matrix_calloc(nTerms, nTerms);   
	double chisq;  
	  
	gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(length, nTerms);   
  
	int res = gsl_multifit_linear (x,  
								   yx,  
								   cx,  
								   cov,  
								   &chisq,  
								   work);  
	  
	int res2 = gsl_multifit_linear (x,  
								   yy,  
								   cy,  
								   cov,  
								   &chisq,  
								   work);  
	  
	printf("-------------------------------------------- \n");  
	  
	  
	double * xptr = gsl_vector_ptr(cx,0);  
	double * yptr = gsl_vector_ptr(cy,0);  
  
	for (int i = 0; i < nTerms; i++){  
		printf("cx %i = %f \n", i, xptr[i]);  
		cxfit[i] =  xptr[i];  
	}  
	  
	for (int i = 0; i < nTerms; i++){  
		printf("cy %i = %f \n", i, yptr[i]);  
		cyfit[i] =  yptr[i];  
	}  
	  
	  
	bBeenFit = true;  
	  
	  
	  
	printf("-------------------------------------------- \n");  
	  
	  
	//return ;  
	  
	  
}  
  

// how to use the remapping:

  
  
float remappedx = (x) * cxfit[0] + (x*x) * cxfit[1] + (y) * cxfit[2] + (y*y) * cxfit[3]   
						+ (x*y) * cxfit[4] + (x*x*x) * cxfit[5] + (y*y*y) * cxfit[6] + cxfit[7];  
  
float remappedy = (x) * cyfit[0] + (x*x) * cyfit[1] + (y) * cyfit[2] + (y*y) * cyfit[3]   
						+ (x*y) * cyfit[4] + (x*x*x) * cyfit[5] + (y*y*y) * cyfit[6] + cyfit[7];  
  

btw, it was the amazing Tim Hanson who sat with me for an hour to explain this and get it to work, who is part of this great project openprosthetics.org – please check it out, an effort to help support research and development of open source prosthetic limbs.

I misread the title and got super excited.

I thought it said “linear regression example w/ GLSL” :slight_smile: Still, nice tip.

Also, OpenCV lets you do this too with CvStatModel.