How to fit Plane (z=ax+by+c) to 3D point cloud data?

Hi, I am trying to do plane fit to 3D point data. Point cloud file is attached. Here is my code I tried using least square method

clc;
clear all;
close all;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%reading %%%%%%%%%
arr = step_mask('Step_scan01_ex.xlsx','Sheet1', 'A:C');
 %subplot(1,3,1)
x=arr(:,1);
y=arr(:,2);
z=arr(:,3)
plot3(arr(:,1),arr(:,2),arr(:,3),'.'); grid on
xlabel('x(mm)'); ylabel('y(mm)'); zlabel('z(mm)');
title('Masked plot');
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%least squares method %%%%%%%%%%%

   [xnum,ynum]=size(arr);
   % 

   for i = 1: xnum 

   xz (i) = x (i) * z (i);
   yz (i) = y (i) * z (i); 
   xy (i) = x (i) * y (i); 
   end 

   Sxz = sum (sum (xz)); 
   Syz = sum (sum (yz)); 
   Sz = sum (sum (z)); 
   Sxy = sum (sum (xy)); 
   Sx =  ynum*sum (x); 
   Sy =  sum (y); 
   Sx2 = sum (x.^2); 
   Sy2 =  sum (y.^2); 

   A = [Sx2 Sxy Sx; Sxy Sy2 Sy; Sx Sy xnum * ynum]; 
   Z = [Sxz; Syz; Sz]; 

   W = A \ Z; 

   w1 = W (1); 
   w2 = W (2); 
   w3 = W (3); 

   l = - w1/(w1^2 + w2^2 + 1)^0.5; 
   m = - w2/(w1^2 + w2^2 + 1)^0.5; 
   n = 1/(w1^2 + w2^2 + 1)^0.5; 
   p = w3/(w1^2 + w2^2 + 1)^0.5;
   %%%%%%%%%%%%%%%%Generating and displaying the obtained plane %%%%%%%%%
   z_2=zeros(7340,1);
   for i = 1: xnum 

    z_2(i)=(-l/n)*x(i)+(-m/n)*y(i)+(p/n);
    i=i+1;
   end 

a=-l/n;
b=-m/n;
c=p/n;
averagev = sum(sum(z_2))/(xnum * ynum); 

   figure;
   plot3(x,y,z_2,'.'); grid on
   xlabel('x'); ylabel('y'); zlabel('z(mm)');
   title('Fitted Plane');

   sa = abs(z-z_2);
   figure;
  plot3(x,y,sa,'.'); grid on
   xlabel('x'); ylabel('y'); zlabel('z(mm)');
   title('Substracted Plane');

Result of this code is not looking fine to me. Could anyone please suggest what is the best way to fit plane with 3D pint data?

Answer 1

[Star Strider]

Try this:

arr = xlsread('Step_scan01_ex.xls');

x=arr(:,1);
y=arr(:,2);
z=arr(:,3);

DM = [x, y, ones(size(z))];                             % Design Matrix
B = DM\z;                                               % Estimate Parameters

[X,Y] = meshgrid(linspace(min(x),max(x),50), linspace(min(y),max(y),50));
Z = B(1)*X + B(2)*Y + B(3)*ones(size(X));

figure(1)
plot3(arr(:,1),arr(:,2),arr(:,3),'.')
hold on
meshc(X, Y, Z)
hold off
grid on
xlabel('x(mm)'); ylabel('y(mm)'); zlabel('z(mm)');
title('Masked plot');
grid on
text(-20, 50, 450, sprintf('Z = %.3f\\cdotX %+.3f\\cdotY %+3.0f', B))

[Swati Jain]

Thank You!

[Star Strider]

As always, my pleasure!

[Manuel Fuelling]

Do i understand the code correct?: Having a fitted plane with the form z=ax+by+d i should have the values for a,b and d in B(1), B(2) and B(3) respectively.

[Star Strider]

Yes.

[Stephen]

Is this fitting mechanism using orthogonal distance regression, instead of just minimizing the error in Z direction?