You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
3D line of best fit
9 views (last 30 days)
Show older comments
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
Accepted Answer
Matt J
on 7 Jun 2013
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);
31 Comments
kanai
on 7 Jun 2013
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
on 7 Jun 2013
Edited: Matt J
on 7 Jun 2013
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
on 22 Jun 2013
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
on 23 Jun 2013
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.
Matt J
on 23 Jun 2013
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
on 23 Jun 2013
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
on 23 Jun 2013
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
on 23 Jun 2013
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
on 23 Jun 2013
Edited: Matt J
on 23 Jun 2013
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
on 24 Jun 2013
Edited: Matt J
on 24 Jun 2013
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
on 24 Jun 2013
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
on 24 Jun 2013
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
on 25 Jun 2013
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
on 27 Jun 2013
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
on 27 Jun 2013
Edited: Matt J
on 27 Jun 2013
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
on 4 Jul 2013
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
Matt J
on 16 Jul 2018
Edited: Matt J
on 16 Jul 2018
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
on 18 Apr 2019
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
on 18 Apr 2019
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
![untitled.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/214785/untitled.png)
Songqiu
on 2 Oct 2020
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
Songqiu
on 2 Oct 2020
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
More Answers (0)
See Also
Categories
Find more on Line Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)