Pointing to my previous email, I've just wanted to let you know the algorithm that I used for the pyramid using Matlab code for computing the apex. The parameters of (a,b,c) for each plane I calculated separately by another algorithm. But the way for calculating stdv. for apex was wrong. YOU solved my problem by your nice algorithm

clear all
clc
C1=load('D:\New Folder\pos114.10.11\p11x1py4\c1p1x1.txt'); % point of plane 1
X=C1(:,1);
Y=C1(:,2);
Z=C1(:,3);
figure(1)
plot3(X,Y,Z,('k.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid 
 
x = C1(:,1);
y = C1(:,2);
a=-0.21131;b=0.87480;c=-0.43598;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(2)
L=plot3(x,y,z,'ko'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','k') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d1=z-C1(:,3);
sdvsurf1=std(zz);
sdv1=std(d1)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
D1=load('D:\New Folder\pos114.10.11\p11x1py4\d1p1x1.txt'); % points of plane 2
X=D1(:,1);
Y=D1(:,2);
Z=D1(:,3);
figure(3)
plot3(X,Y,Z,('r.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = D1(:,1);
y = D1(:,2);
a=0.280048226488374;b=0.959535677;c=0.02939857;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(4)
L=plot3(x,y,z,'ro'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
%set(L,'Markerfacecolor','r') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d2=z-D1(:,3);
sdvsur2=std(zz);
sdv2=std(d2)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
F1=load('D:\New Folder\pos114.10.11\p11x1py4\f1p1x1.txt'); % points of plane 3
X=F1(:,1);
Y=F1(:,2);
Z=F1(:,3);
figure(7)
plot3(X,Y,Z,('g.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = F1(:,1);
y = F1(:,2);
a=-0.15202644;b=0.830514033;c=0.535849234;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(8)
L=plot3(x,y,z,'g+'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
%set(L,'Markerfacecolor','g') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d4=z-F1(:,3);
sdvsur4=std(zz);
sdv4=std(d4)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
E1=load('D:\New Folder\pos114.10.11\p11x1py4\e1p1x1.txt'); %point of plane 4
X=E1(:,1);
Y=E1(:,2);
Z=E1(:,3);
figure(5)
plot3(X,Y,Z,('b.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = E1(:,1);
y = E1(:,2);
a=-0.626881329;b=0.773774228;c=0.09106725;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(6)
L=plot3(x,y,z,'b*'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','b') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid (min(X):0.01:max(X),min(Y):0.01:max(Y),min(Z):0.01:max(Z)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d3=z-E1(:,3);
sdvsur3=std(zz);
sdv3=std(d3)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
% By rotating the surface, you can see that the points lie on the plane
% Also, if we multiply both sides of the equation in the title by 4,
% we get the equation in the comment on the third line of this example
 
sdvapex=sqrt(sdv1^2+sdv2^2+sdv3^2+sdv4)/4
 
C1=load('D:\New Folder\pos114.10.11\p11x1py4\c1p1x1.txt');
X=C1(:,1);
Y=C1(:,2);
Z=C1(:,3);
figure(9)
plot3(X,Y,Z,('k.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid 
hold on
 
D1=load('D:\New Folder\pos114.10.11\p11x1py4\d1p1x1.txt');
X=D1(:,1);
Y=D1(:,2);
Z=D1(:,3);
plot3(X,Y,Z,('ro'))
 
E1=load('D:\New Folder\pos114.10.11\p11x1py4\e1p1x1.txt');
X=E1(:,1);
Y=E1(:,2);
Z=E1(:,3);
plot3(X,Y,Z,('b*'))
 
F1=load('D:\New Folder\pos114.10.11\p11x1py4\f1p1x1.txt');
X=F1(:,1);
Y=F1(:,2);
Z=F1(:,3);
plot3(X,Y,Z,('g+'))