3D line of best fit

I have about 50000 points with x,y,z data spread in 3 coloumns in excel. I have managed to create a plane of best fit. (by creating a comma delimited file, importing it, setting each column as a variable, then using the SFTOOL function. This gives me a plane of best fit through the data when I select the plane to be polynomial with degrees x=1 and y=1).

What I cannot seem to do is create a line of best fit through this data. This line must be a straight line and I would like to know the parametric equation of this line.

I would be most grateful if someone could advise me (a complete matlab beginner) how to go about this.

Many thanks

regards

kanai

Answer 1

[Matt J]

Suppose the rows of xyz are your data points,

 r0=mean(xyz);
 xyz=bsxfun(@minus,xyz,r0);

 [~,~,V]=svd(xyz,0);

The parametric equation is

 r(t) = r0 + t*V(:,1);

[kanai]

Many thanks for getting back to me Matt. I really appreciate it.

At the risk of sounding like an idiot: When I put r0=mean(xyz) into the command window it comes up with an error (??? Undefined function or variable 'xyz'.)

I think I have specified my x y z data by right clicking the columns in the variable editor and selecting create a variable. Then choosing X or Y or Z.

Is there something I am doing wrong?

many thanks

[Matt J]

Suppose the rows of xyz are your data points,

To clarify this part, I meant that you should import your data into an Nx3 matrix called xyz where each row of xyz is one of your 3D data points and N is 50000 or however many points you said you had.

[kanai]

Sorry Matt

I didn't get an email informing me you had responded or I would have gotten back to you far sooner.

Many thanks for getting back to me. I will give you some background to myself and my project.

I am a trainee surgeon in the UK. My 9-5 job is working in the hospital. I am also pursuing a PhD on the side. Unfortunately this makes me a complete computer idiot.

I have managed to play around with matlab to sort out one of my problems (namely a plane of best fit)

I would be very grateful if could give me a step-by-step guide on how to go from xyz coordinates in 3 separate columns in an excel file (comma/tab delimited) to a line of best fit.

I apologise for not being a computer whizz! I also notice your interests are medical image processing!

yours sincerely

kanai

[Matt J]

I'm a bit confused then. How, as in your initial description, did you manage to import from a comma delimited file and use sftool if you are "a complete computer idiot"?

Just import the coordinates again as you did before. Or, if they are in Excel form now (how did they get there from the comma delimited file?), then import them using XLSREAD.

[kanai]

Alot of trial and error!! I have no problem uploading data, but how do you upload to an Nx3 matrix?

Does the data imported need to be set out in rows rather than columns for the matrix?

[Matt J]

xyz=xlsread(filename.xls,range)

will read from the 'range' of cells in the spreadsheet into a matrix xyz. The matrix xyz will have the same shape as in the spreadsheet. If the data is Nx3 in the spreadsheet, you will have what you need to run my earlier code. If the data is 3xN in the spreadsheet, xyz will also read in as 3xN, but you can simply transpose it:

xyz=xyz.'

[kanai]

I have uploaded the comma delimited file. I have opened the file with 3 columns. I have created a variable from each column by right clicking on the top box and creating a variable. I have named these xyz. I have tried all the commands you mention but i keep getting undefined function or variable xyz. When i type x or y or z and press enter, i get my data.

What am I doing wrong?

[Matt J]

I have named these xyz

No, I don't think you did. I think what you really mean is that you have named the different variables x, y, and z.

If so, then make them into a single Nx3 array as follows

 xyz=[x(:), y(:), z(:)]

Now, run my code.

[kanai]

I think it is really close now and I thank you for all the help

I have done:

1) xyz=[x(:), y(:), z(:)]

2) r0=mean(xyz);

3) xyz=bsxfun(@minus,xyz,r0);

4) [~,~,V]=svd(xyz,0);

Then when I do r(t) = r0 + t*V(:,1); I get ??? Undefined function or variable 't'. When I do [~,~,V]=svd(xyz,0) I get

V= Columns 1 through 2

    0.4475   -0.1748
   -0.4678    0.7901
    0.7622    0.5876

Column 3

    0.8771
    0.3962
   -0.2717

How do I convert this to for example Line x= 2+3t y=3+2t z=4+3t etc.

Many thanks

[Matt J]

Then when I do r(t) = r0 + t*V(:,1); I get ??? Undefined function or variable 't'.

That wasn't meant to be a line of code. Line 4) was the final line of code.

You're done now. You have your parametric equation. r0 is the reference point on the line and V(:,1) (the first column of V) is the direction vector.

Line x= 2+3t y=3+2t z=4+3t etc.

So, in other words

    x=r0(1)+t*V(1,1)
    y=r0(2)+t*V(2,1)
    z=r0(3)+t*V(3,1)

[kanai]

Just to clarify

So if

r0 =

    Columns 1 through 2

      0.9339  117.6053

    Column 3

    -90.2660

and

V =

      0.4475   -0.1750    0.8770
     -0.4679    0.7899    0.3964
      0.7621    0.5877   -0.2716

Then my parametric equation is

    x= 0.9339 + 0.4475t

    y= 117.6053 - 0.4679t

    z= -90.266 + 0.7621t

Is this correct? If so what is the V values for column 2 and 3 for?

I really can't thank you enough mate! Once this is clarified then I will consider the question answered.

[Matt J]

Is this correct? If so what is the V values for column 2 and 3 for?

Yes, that's correct. The other columns of V are not useful for your current task, but there is no way to tell svd() to give you just the one column you're interested in.

[kanai]

I have accepted the answer. Thanks a million.

Is there a way of plotting all of my points and displaying this line through the data set as well.. I can set up another question if you would prefer to answer it that way?

[Matt J]

Yes, you can use plot3()

[Matt J]

and scatter3()

[kanai]

Thanks. I'd be most appreciative if you provide me with the code I need to type in to get the 3d line and scatter chart on the same graph

Is there something I need to type in the brackets?

[kanai]

I've just had a play and

scatter3 (x,y,z)

or

plot3 (x,y,z)

creates the scatter/plot of my points. Its now just the matter of putting the line through the points that I am unable to do

[Matt J]

Here is a simplifed example

%Fake data
xdata=(0:10) + rand(1,11);
ydata=(20:30) + rand(1,11);
zdata=(30:40) + rand(1,11);

t=0:10;
xfit=0+t;
yfit=20+t;
zfit=30+t;

hold on
 scatter3(xdata,ydata,zdata,'b*')
 plot3(xfit,yfit,zfit,'r')
hold off

[kanai]

thanks, just got back from a LONG shift. Will give it a go tomorrow and let you know how I got on

[kanai]

I have got to where we were last time with V values and r0 values.

do i need to put

xdata=(0:10) + rand(1,11);

ydata=(20:30) + rand(1,11);

zdata=(30:40) + rand(1,11);

into my code or is that just an example.

2) Should I rename my variables xdata, ydata and zdata for this to work?

3) Where you have put xfit/yfit/zfit. Do i put the parametric equations for my line?

4) What is 'b*'?

The main confusion I have is the purpose of

xdata=(0:10) + rand(1,11); ydata=(20:30) + rand(1,11); zdata=(30:40) + rand(1,11);

[Matt J]

do i need to put ... do i need to put

It's just an example, which is why I labeled it "Fake data"

2) Should I rename my variables xdata, ydata and zdata for this to work?

You could do that or you could rewrite the call to scatter3 with whatever variable names you are using now.

3) Where you have put xfit/yfit/zfit. Do i put the parametric equations for my line?

Yes.

4) What is 'b'?

It makes the scatter plot mark the plotted points with blue stars.

I have a suspicion that you are not aware of the DOC command. If you were to run "doc scatter3", for example, you will be shown documentation about how these commands work, with examples and pictures. It would tell you things like what the 'b*' is doing in the code and, more generally, help you adapt examples to your needs...

[kanai]

Whenever I type my

xfit=0.9339+0.4475t i get

"Error: Unexpected MATLAB expression."

I have defined t above this with t=0:10;

Even if I start a new window and try just

t=0:10;

xfit=0.9339 + 0.4475t

I get the error. What am I doing wrong?

many thanks

[kanai]

I figured it out... you have to put 0.4475*t

[duong Mai]

Hi, Matt J. Could you tell me what is the name of your fitting method?

[Matt J]

Hi duong, I call it "Steve's method" (my friend Steve showed it to me). However, it might really be a form of total least squares or principal component analysis.

[P_L]

Hi Matt,

I have looked at using your method to fit a line tyo my 3-D data.

Please see my code:

lat_3d= [30.3894980000000;30.3824250000000;30.3886690000000;30.3878380000000;30.3953250000000;30.3874200000000;30.3807640000000;30.3924160000000;30.3874200000000;30.3811820000000;30.3919950000000;30.3928260000000;30.3903320000000;30.3836720000000;30.3774440000000];%new_3d_model_results(:,2);
lon_3d= [40.7508620000000;40.7486390000000;40.7489620000000;40.7476950000000;40.7515010000000;40.7486440000000;40.7461050000000;40.7486490000000;40.7486440000000;40.7451550000000;40.7514980000000;40.7521320000000;40.7502300000000;40.7495900000000;40.7410360000000];%new_3d_model_results(:,1);
depth_3d=[10.0097660000000;8.04296900000000;8.97460900000000;7.98085900000000;9.51289100000000;8.04296900000000;10.1132810000000;10.0304690000000;7.96015600000000;7.98085900000000;9.47148400000000;7.77382800000000;8.76757800000000;9.34726600000000;7.79453100000000]; %new_3d_model_results(:,3);

%% plot 3d 
scale = 111; % km/deg
figure 
%plot
light_p=[204 153 255]/255;
% h=[30.3 30.45 40.7 40.8];
% lat_lon_proportions(h) % refernced at top of script 

% plot events 
hold on 
p10=plot3(lat_3d,lon_3d,depth_3d/scale,'o','MarkerFaceColor',light_p,'MarkerEdgeColor', 'k','MarkerSize',12); %QPS changed/ Initial Velocity Model
hold on 
grid on 
%axis equal
ax = gca;
ax.ZDir = 'reverse';
zVal = str2double(ax.ZTickLabel)*scale;
ax.ZTickLabel = num2str(zVal);
xlabel('Longitude (deg, E)');
ylabel('Latitude (deg, N)')
zlabel('Depth (km)')
zval = [2:2:20]';
ax.ZTick = zval/scale;
ax.ZTickLabel = num2str(zval);
hold on 

set(gca,'FontName','Helvetica','FontSize',20);
hold on 
xlim([30.3 30.45])
ylim([40.7 40.8])
hold on 
axis([30.3 30.45 40.7 40.8]);
view(-37.5, 4)


% import data into matrix 
    xyz=[lon_3d(:) lat_3d(:) depth_3d(:)];
    r0=mean(xyz);
    xyz=bsxfun(@minus,xyz,r0);
    [~,~,V]=svd(xyz,0);
    
    t=1:3
    r(t) = r0 + t*V(:,1);
    hold on 
    xfit=r0(1)+t*V(1,1)
    yfit=r0(2)+t*V(2,1)
    zfit=r0(3)+t*V(3,1)
    hold on 
    plot3(xfit,yfit,zfit,'r')
    

when I/ you run it, I can't get the line to plot correctly - please could you help me correct this?

Many thanks

[Matt J]

    xyz=[lon_3d(:) lat_3d(:) depth_3d(:)/scale];
    r0=mean(xyz);
    [~,~,V]=svd(bsxfun(@minus,xyz,r0),0);

    t=linspace(-.1/3,.1/3,10);
    xfit=r0(1)+t*V(1,1);
    yfit=r0(2)+t*V(2,1);
    zfit=r0(3)+t*V(3,1);

p10=plot3(xyz(:,1),xyz(:,2),xyz(:,3),'o',...
      'MarkerFaceColor',light_p,'MarkerEdgeColor', 'k','MarkerSize',12); %QPS changed/ Initial Velocity Model
 
 line(xfit,yfit,zfit)    
 grid on
 axis vis3d

[Songqiu]

Hi Matt,

How can we get the xyz values of the projected location to the fitted line for each point? (i.e. the point on the fitted line which has the shortest distance to the corresponding point).

Thanks a lot.

Songqiu

[Matt J]

This may help,

https://gamedev.stackexchange.com/questions/72528/how-can-i-project-a-3d-point-onto-a-3d-line

[Songqiu]

Hi, Matt,

Thank you very much for your reply.

I have written the following code for calculating the projections of multiple points using a for loop.

Could you please tell me how to avoid using the for loop?

Many thanks.

a = [0,  0, 0];
b = [1, 1, 1];
p = [1, 1, 0];
for i = 2 : 50
    p = cat(1, p, [i i 0]);
end

pp = p;
xyz = zeros(size(pp, 1), 3);
for i = 1 : size(pp, 1)
    p = pp(i, :);
    ap = p-a;
    ab = b-a;
    xyz(i, :) = a + dot(ap, ab)/dot(ab, ab) * ab;
end

[Matt J]

ab=(b-a)/norm(b-a);

xyz=(pp-a)*ab(:)*ab + a;