point cloud with matlab

Hello Matlab,

I really want to ask regarding 3D Point cloud. we have points with X, Y, Z. We are totally new at this work. With MatLab we want to do 2D point cloud.

Any one have reference or solution etc?

[Stephan]

As you can see by following the discussion, it would be helpful to have some more informations or maybe a plot of the data in 3D.

[Bruno Luong]

What is best-fit-algorithm? I'm not familiar with it (why people pretentiously call something generally "best" is beyond my understanding).

If you know the rotational axis (by 2 points or whatever), a simple solid-geometry transformation can match it to z-axis. This is independent question.

And why the fit method has something to do with the coordinate transformation?

I don't get why all this questions?

[Image Analyst]

"are almost similar to solution of this task"??? You mean there is an official solution of this somewhere? I don't remember you telling us this. What is the "correct" solution of this task? Can you share it with us?

[Image Analyst]

OK, but I agree with Bruno that we don't know what the "Best Fit" algorithm is, if there even is an algorithm by that name. The "least squares method" is the "maximum likelihood" solution so one might call that the, or a, best fit algorithm. I think least squares is "best" only assuming the residuals have a normal distribution, which is not always the case. Of course there are other algorithms that follow different approaches depending on what you know or want to assume (like noise distribution).

Answer 1

[Bruno Luong]

OP's question

"I really want to ask regarding 3D Point cloud. we have points with X, Y, Z. We are totally new at this work. With MatLab we want to do 2D point cloud.

Any one have reference or solution etc?"

xyz=xlsread('PointClouds.xlsx');

r=sqrt(sum(xyz(:,1:2).^2,2));
z=xyz(:,3);

P = polyfit(z,r,5);
zi = linspace(min(z),max(z),50);
ri = polyval(P,zi);
theta = linspace(0,2*pi,100);
[THETA,Z] = meshgrid(theta,zi);
R = polyval(P,Z);
[X,Y,Z] = pol2cart(THETA,R,Z);

close all
figure();
subplot(2,1,1);
h = plot(r,z,'.b',-r,z,'.b',ri,zi,'r','linewidth',2);
axis([-80 80 -10 100])
axis equal
legend(h([1 3]),'reduction 2D','fit')
subplot(2,1,2);
surf(X,Y,Z);
axis equal

[Stephan]

Very nice! +1

[Bruno Luong]

Attached are data that have been removed by OP

Answer 2

[Image Analyst]

Try this:

% Read data from file.
data = xlsread('PointClouds.xlsx')
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
% Plot 3-D point cloud.
subplot(1, 2, 1);
scatter3(x, y, z, 1);
grid on;
xlabel('X', 'FontSize', 15);
ylabel('Y', 'FontSize', 15);
zlabel('Z', 'FontSize', 15);
title('3-D Point Cloud', 'FontSize', 15);

% Plot 2-D point cloud.  Ignore Z.
subplot(1, 2, 2);
scatter(x, y, 4, 'filled');
grid on;
xlabel('X', 'FontSize', 15);
ylabel('Y', 'FontSize', 15);
title('2-D Scatterplot', 'FontSize', 15);

[M.S. Khan]

Hi Matlab community,

i am using this code for getting two overlaped or intersected spheres. Can anyone guide me how to use above code to get point cloud of such figure. Thanks in advance for all cooperation and guidance.

[X, Y, Z]=sphere;
a = [1 1 1 1.2; 2.293,1,1,1.8];
%Draw sphere #1
s1=surf(X*a(1,4)+a(1,1), Y*a(1,4)+a(1,2), Z*a(1,4)+a(1,3),'FaceColor', [1 1 1],'edgecolor','none','FaceAlpha',0.6);
light
axis equal
set(gca,'Color','w')
hold on
% Draw sphere # 2
s2=surf(X*a(2,4)+a(2,1), Y*a(2,4)+a(2,2), Z*a(2,4)+a(2,3),'FaceColor', [0 1 0],'edgecolor','none','FaceAlpha',0.5);
daspect([1 1 1])
view(30,20)
xlabel('x')
ylabel('y')
zlabel('z')
x1range1 = -0.2;
x1range2 = 4.093;
y1range1 = -0.8;
y1range2 = 2.8;
z1range1 = -0.8;
z1range2 = 2.8;

Answer 3

[Stephan]

Hi,

i would expect this to work:

[~,x_2D,y_2D] = cart2pol(x,y,z);
scatter(x_2D,y_2D)

The result should be that all the points from 3D space are converted into polar coordinates, the angle does not matter. It is thus pretended that all points lie in one plane. For the rotational symmetry you specify with respect to z, this should lead to the correct result:

% Read data from file.
data = xlsread('PointClouds.xlsx');
x_org = data(:, 1);
y_org = data(:, 2);
z_org = data(:, 3);

% Plot 2-D point cloud using cart2pol function.  
[~,x_2D,z_2D] = cart2pol(x_org,y_org,z_org);
scatter(x_2D,z_2D,'+')

% Fit a function to the 2D points an plot
fun = @(x,x_2D)(x(1).*x_2D.^4 + x(2).*x_2D.^3 + x(3).*x_2D.^2 + x(4).*x_2D + x(5))./(x_2D.^2 + x(6).*x_2D + x(7));
x = lsqcurvefit(fun,[3 21 13 -9 22 -1 0.7],x_2D,z_2D);
z_fit = fun(x,x_2D);
hold on
scatter(x_2D,z_fit,15,'or','filled')
legend({'2D Data', 'Fit'})
hold off

results in:

Note that this fit is a rational function, which i found that it fits best. Disadvantage: Optimization Toolbox is needed for this code. The fit function is:

f(x) = (a*x^4 + b*x^3 + c*x^2 + d*x + e)/(x^2 + f*x + g)

with results for a-g:

   1.0e+04 *

   0.000003831224870
  -0.000477260475870
   0.018345571140362
  -0.159439995510358
  -1.991256255238212
  -0.010298052678467
   0.268478271805655

Best regards

Stephan

[Image Analyst]

Clever. How did you come up with that function, rather than just a simple 5th order polynomial like Bruno used? I don't have that toolbox - is it one of the options in that toolbox?

What is the residual for that function compared to the 5th order polynomial?

[Bruno Luong]

Here is my fit. I use the polynomial but r=P(z) and not the opposite.

Let you judge which one is better (look at the begining)

[Stephan]

Yes it is an option of curve fitting toolbox see here. I just wanted to show another approach that might be useful. The results for the residuals of the different approaches will tell us the truth. Until now i did not compare the results.

I think we all learn from each other by sharing different approaches and this is resulting in a wide range of useful answers. If i look at Brunos graph, i think his one seems to be a bit better fit at the start.

Here are the results of:

[fitobject,gof] = fit(x_2D,z_2D,'rat42')

fitobject = 

     General model Rat42:
     fitobject(x) = (p1*x^4 + p2*x^3 + p3*x^2 + p4*x + p5) /
               (x^2 + q1*x + q2)
     Coefficients (with 95% confidence bounds):
       p1 =     0.03885  (0.03776, 0.03994)
       p2 =      -4.857  (-5.048, -4.666)
       p3 =       188.5  (176.3, 200.7)
       p4 =       -1728  (-2065, -1390)
       p5 =  -1.859e+04  (-2.2e+04, -1.517e+04)
       q1 =      -102.9  (-103, -102.8)
       q2 =        2681  (2674, 2687)

gof = 

  struct with fields:

           sse: 2.0329e+04
       rsquare: 0.9980
           dfe: 9993
    adjrsquare: 0.9980
          rmse: 1.4263

Compared to the same with 5th order polynomial:

[fitobject,gof] = fit(x_2D,z_2D,'poly5')
         
fitobject = 

     Linear model Poly5:
     fitobject(x) = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6
     Coefficients (with 95% confidence bounds):
       p1 =  -2.127e-05  (-2.169e-05, -2.085e-05)
       p2 =    0.004077  (0.003989, 0.004164)
       p3 =     -0.3018  (-0.3089, -0.2947)
       p4 =       10.85  (10.56, 11.13)
       p5 =      -188.5  (-193.9, -183)
       p6 =        1257  (1216, 1298)

gof = 

  struct with fields:

           sse: 4.2389e+04
       rsquare: 0.9958
           dfe: 9994
    adjrsquare: 0.9958
          rmse: 2.0595

It seems to be better using rat42 - but Bruno did it this way:

[fitobject,gof] = fit(z_2D,x_2D,'poly5')

itobject = 

     Linear model Poly5:
     fitobject(x) = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6
     Coefficients (with 95% confidence bounds):
       p1 =   2.627e-10  (-2.957e-11, 5.551e-10)
       p2 =  -5.946e-08  (-1.258e-07, 6.852e-09)
       p3 =   7.106e-05  (6.59e-05, 7.623e-05)
       p4 =    -0.01178  (-0.01194, -0.01162)
       p5 =      0.8003  (0.7986, 0.802)
       p6 =       30.51  (30.5, 30.52)

gof = 

  struct with fields:

           sse: 832.7220
       rsquare: 0.9990
           dfe: 9994
    adjrsquare: 0.9990
          rmse: 0.2887
          

This is a much better result. Now try the same with rat54 option (rat42 performs really bad in this case):

[fitobject,gof] = fit(z_2D,x_2D,'rat54')

fitobject = 

     General model Rat54:
     fitobject(x) = 
               (p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6) /
               (x^4 + q1*x^3 + q2*x^2 + q3*x + q4)
     Coefficients (with 95% confidence bounds):
       p1 =       101.5  (-589.1, 792.1)
       p2 =  -1.781e+04  (-1.394e+05, 1.038e+05)
       p3 =    1.25e+06  (-7.289e+06, 9.789e+06)
       p4 =   4.475e+07  (-2.605e+08, 3.5e+08)
       p5 =  -6.648e+07  (-5.2e+08, 3.87e+08)
       p6 =   2.388e+07  (-1.39e+08, 1.868e+08)
       q1 =       176.9  (-1388, 1742)
       q2 =   1.524e+06  (-8.874e+06, 1.192e+07)
       q3 =  -2.199e+06  (-1.72e+07, 1.28e+07)
       q4 =   7.827e+05  (-4.556e+06, 6.122e+06)

gof = 

  struct with fields:

           sse: 833.7492
       rsquare: 0.9990
           dfe: 9990
    adjrsquare: 0.9990
          rmse: 0.2889

We can see that Brunos way performs still a little better and additionally with a much easier approach!

Thank you Bruno for this knowledge-profit!

[Bruno Luong]

Nice analysis

Answer 4

[Image Analyst]

Why not simply ignore the z value? It would be just like you're looking down along the Z axis so that all points are like mapped to the x-y plane?

If you want it to be symmetric about the z axis (out of the plane), simply subtract the mean x from all the x values, and the mean y from all the y values. Then the cluster will be centered at the origin (where x and y axes cross) instead of off away somewhere.

[Stephan]

The way i understand this question is that the points in cartesian coordinates describe a rotational symmetric curve with respect to the z-axis. I imagine a bottle of wine (mmh... ;-).

So if i am right, the view from the top to x-y-plane would give a lot of circles with different diameters. I dont see how to fit a curve with this information.

Maybe im wrong, but this is how i understand the question.

[Image Analyst]

Not what I envisioned. I didn't see anything about the x,y,z being coordinates on one or more curves. I'm envisioning a "3D point cloud" being basically like a shotgun blast of unrelated points. Posting a screenshot of the 3-D scatterplot would certainly help.

[Stephan]

?