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

[John Pye]

G'day to you too...

Thanks for the reply, John. I think you're right. In general a plane
requires the 'd' parameter, but only really for the degenerate case
where the plane passes through the origin. For any other case, the
equation can be rearranged such that the 'd' value is replaced by a '1'
on the RHS.

This lead me to wonder how I can implement a plane-fitting routine that
is robust and general. I think that the correct approach might be to
offset all the points by the centroid, and then fit the equation
ax+by+cz=0. There is a theorem about the centroid always lying on the plane.

So now I tried to implement this using GSL but it seems that GSL lacks a
function like 'gsl_multifit_mul'. When I try to fit my plane equation
using gsl_multifit_linear, with y[i]=0 for all i, I get the predictable
result that a=0, b=0, c=0. So I don't think GSL can do what I need it to
do. Or perhaps there is a cunning workaround?

Thoughts?

Cheers
JP

John Gehman wrote:
> G'day John,
>
> I haven't looked at the examples all that closely, but just quickly, I
> believe a general equation for a plane is ax+by+cz+d=0, i.e. in your
> equation you're effectually fixing d at -1, while I think it should be
> d=0 (zero) in your example?
>
> Cheers,
> john
>
> ---------------------------------------------------------
>
>
> John Gehman
>
> Research Fellow
>
> School of Chemistry
>
> University of Melbourne
>
> Australia
>
>
> On 20/11/2007, at 12:53 PM, John Pye wrote:
>
>> Hi all
>>
>> I have a question about the use of the gsl_multifit_linear routine that
>> perhaps is a question more about geometry/algebra than coding, but I'm
>> not sure.
>>
>> I want to construct a routine that fits a plane (in three dimensions)
>> through a set of data points (x,y,z). I have set up gsl_multifit_linear
>> to fit the plane equation a*x + b*y + c*z = 1to my data, and for most
>> cases that seems to work OK. However there are a few degenerate cases
>> that don't work, and I'm trying to work out what I should do. Is there a
>> better equation that describes a plane?
>>
>> The following example shows what happens when I try to fit a plane to
>> the points (0,0,0), (1,0,0), (1,1,1), (0,1,1). These points represent an
>> inclined rectangle with one side on the positive x-axis.
>>
>> I would expect the fit to these points to be 0*x + INF*y - INF*z = 1, or
>> else for it to give an error. But instead, gsl_multifit_linear gives me
>> the result 0.666*x + 0.333*y + 0.333*z = 1, which is wrong.
>>
>> I am obviously making a logic error here, but I'm a bit stuck on what it
>> is, and what the correct general way of working around it would be? Can
>> anyone here offer any suggestions?
>>
>> Cheers
>> JP
>>
>>
>> // for the fitting of a plane through a group of points, testfit.cpp
>> #include <gsl/gsl_vector.h>
>> #include <gsl/gsl_matrix.h>
>> #include <gsl/gsl_multifit.h>
>> #include <iostream>
>>
>> using namespace std;
>>
>> int main(void){
>>
>>     int n = 4;
>>     int num_params = 3;
>>
>>     gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,
>> num_params);
>>
>>     gsl_vector *y = gsl_vector_calloc(n);
>>     for(unsigned i=0; i < n; ++i){
>>         gsl_vector_set(y,i,1.0);
>>     }
>>
>>     gsl_matrix *X = gsl_matrix_calloc(n, num_params);
>>     gsl_vector *c = gsl_vector_calloc(num_params);
>>     gsl_matrix *cov = gsl_matrix_calloc(num_params, num_params);
>>
>> #define SET(I,A,B,C) \
>>     gsl_matrix_set(X,I,0,A); gsl_matrix_set(X,I,1,B);
>> gsl_matrix_set(X,I,2,C)
>>
>>     SET(0, 0,0,0);
>>     SET(1, 1,0,0);
>>     SET(2, 1,1,1);
>>     SET(3, 0,1,1);
>>
>>     double chisq;
>>     int res = gsl_multifit_linear(X,y,c,cov,&chisq,work);
>>
>>     cerr << "FIT FOUND:" << endl;
>>     for(int i=0;i<3;++i){
>>         cerr << "  c[" <<i << "]=" <<gsl_vector_get(c,i) << endl;
>>     }
>>
>> }
>>
>>
>> _______________________________________________
>> Help-gsl mailing list
>> address@hidden <mailto:address@hidden>
>> http://lists.gnu.org/mailman/listinfo/help-gsl
>