help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-gsl] puzzled with result from gsl_multifit_linear...


From: John Pye
Subject: Re: [Help-gsl] puzzled with result from gsl_multifit_linear...
Date: Tue, 20 Nov 2007 17:36:36 +1100
User-agent: Thunderbird 2.0.0.6 (X11/20071022)

Hi Patrick

Patrick Alken wrote:
>> Ok, I've used the gsl_linalg_SV_decomp function to obtain the singular
>> value decomposition, and indeed I have a very small (1e-33) value in my
>> vector 'S'.
>>
>> I'm still feeling my way around GSL: what do you suggest I do next to
>> find the correct fitting plane? Can I use the remaining values in 'S' to
>> form my solution? Or do I need to set up another linear fit with a
>> reduced set of variables?
>>     
>
> Well if you find a zero value in position 'i' in S, then the solution
> is given by the column i of the matrix V as output by the svd
> routine.
>
> Some more info:
> http://en.wikipedia.org/wiki/Singular_value_decomposition#Range.2C_null_space_and_rank
>
> So if you find a tiny singular value, just set the coefficient vector
> equal to the corresponding column of V. Otherwise use the normal
> gsl_multifit_linear.
>
> In the meantime we can start thinking if its possible to modify
> gsl_multifit_linear to handle this case.
>   

OK, I managed to implement this using the gsl_linalg_SV_decomp and
gsl_linalg_SV_solve functions as following. In case it's useful to
anyone else. It nicely picks up the case of points being colinear as well.

Cheers
JP


static void centroid(gsl_matrix *X, double n, double num_params, double
*cent){
/* prob there is a more GSL-ish way of doing this... */
int i,j;
for(j=0;j<num_params;++j)cent[j]=0;
for(i=0;i<n;++i){
for(j=0;j<num_params;++j)cent[j]+=gsl_matrix_get(X,i,j);
}
for(j=0;j<num_params;++j)cent[j]/=n;
}

Plane
PlaneFit::leastSquaresPlane(){

unsigned n = G.size();

if(n < 3){
throw runtime_error("Not enough points for PlaneFit");
}

int num_params = 3;

gsl_matrix *X = gsl_matrix_calloc(G.size(), num_params);
gsl_vector *c = gsl_vector_calloc(num_params);

// load data into the matrix 'X'
int ii=0;
for(Group::const_iterator i=G.begin(); i<G.end(); ++i){
Vector xi = (*i)->pos();
gsl_matrix_set(X,ii,0, xi.x);
gsl_matrix_set(X,ii,1, xi.y);
gsl_matrix_set(X,ii,2, xi.z);
++ii;
}

// calculate the centroid
double cent[3];
centroid(X, n, num_params, cent);

// offset the points by the centroid vector
for(int i=0;i<n;++i){
for(int j=0;j<num_params;++j){
double v = gsl_matrix_get(X,i,j);
gsl_matrix_set(X,i,j, v - cent[j]);
}
}

// calulate the fit using SVD directly (so we catch degenerate cases)
gsl_matrix *V = gsl_matrix_calloc(num_params,num_params);
gsl_vector *S = gsl_vector_calloc(num_params);
int res = gsl_linalg_SV_decomp_jacobi(X,V,S);

// check for zero values in the vector 'S', this indicate redundant
terms in the equation
unsigned have_zero=0;
unsigned j_zero;
double tol = 1e-10;
for(int j=0;j<num_params;++j){
if(gsl_vector_get(S,j) < tol){ // guaranteed greater than zero so no 'fabs'
have_zero++;
j_zero = j;
}
}
if(have_zero > 1){
// perhaps this case occurs when the points are colinear
throw runtime_error("PlaneFit: Points do not define a plane (MULTIPLE
ZEROS IN 'S' VECTOR)");
}else if(have_zero==1){
//cerr << "SINGULAR!" << endl;
for(int j=0;j<num_params;++j){
gsl_vector_set(c,j, gsl_matrix_get(V,j,j_zero));
}
}else{
// the general case, use the SVD to solve X*c=0
gsl_vector *y = gsl_vector_calloc(n);
for(int i=0; i<n; ++i){
gsl_vector_set(y,i,0.0);
}
int res = gsl_linalg_SV_solve(X, V, S, y, c);

gsl_vector_free(y);
}

// resulting plane fit is 0 = (x - x_0)*c[0] + (y-y_0)*c[1] + (z-z_0)*c[2]
/*cerr << "FIT FOUND:" << endl;
for(int i=0;i<3;++i){
cerr << " c[" <<i << "]=" <<gsl_vector_get(c,i) << endl;
}*/

// create plane using direction vector 'd' and location 'p' (the centroid)
Vector d(gsl_vector_get(c,0),gsl_vector_get(c,1),gsl_vector_get(c,2));
Vector p(cent[0], cent[1], cent[2]);
Plane P(p,d);
//cerr << "CENTROID = " << p << endl;

//cerr << "FIT PLANE = " << P << endl;

gsl_matrix_free(X);
gsl_vector_free(c);
gsl_matrix_free(V);
gsl_vector_free(S);

return P;
}






reply via email to

[Prev in Thread] Current Thread [Next in Thread]