Hi Meche,

Thanks a lot for all your valuable information and guides. I'm now is about to finish my thesis I have to submit my thesis by January- February 2013. surely, once, I finish my thesis, I will send a copy of it to you.
The problem of the precision of the apex was remained until I contacted with you, and you gave me a good indication of how precisely estimate. I did a program through matlab code using plane equation to determine the apex. The result of apex coordinates is quite similar to your result. But I couldn't determine the precision of the apex. A copy of the second part of the program(computation of the apex) I attached below.

I also would like to let you know for other targets such as cone and sphere I used JAG3D, as another method for fitting cone and sphere. Is this software is also refer to your self? if yes, can you write the source as you did for pyramid. This software I used since February 2012.

Can I ask if possible to send to me the program of JAG3D for fitting cone and sphere instead of software in order to be able to use through matlab. As the same as you did for pyramid.

Many thanks indeed. You are really a good person I saw it yet during the period of my PhD study.

THEINTERSECTIONPOINT =
 
   -8.7482
   10.2641
   -0.3286
 
 clear all
clc
C1=load('D:\New Folder\pos114.10.11\p11x1py4\c1p1x1.txt');
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.0193;b=0.0797;c=-0.0396;
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');
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.0373;b=0.1293;c=0.0040;
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))
 
E1=load('D:\New Folder\pos114.10.11\p11x1py4\e1p1x1.txt');
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.0468;b=0.0576;c=0.0068;
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
 
F1=load('D:\New Folder\pos114.10.11\p11x1py4\f1p1x1.txt');
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.0158;b=0.0857;c=0.0551;
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))
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+'))
b=ones(length(C1),1);
par1=inv(C1'*C1)*C1'*b % first parameters
standreddivationofsurfaceC=std(C1);
%the dark grey coulour, normal tape
 
b=ones(length(D1),1);
par2=inv(D1'*D1)*D1'*b  %second parameters
standreddivationofsurfaceD=std(D1);
 
%the light grey colour, high reflectivity tape
b=ones(length(E1),1);
par3=inv(E1'*E1)*E1'*b  %third parameters
standreddivationofsurfaceE=std(E1);
%the green coulour, normal tape
 
b=ones(length(F1),1);
par4=inv(F1'*F1)*F1'*b  %second parameters
standreddivationofsurfaceF1=std(F1)
Q4=[par1';par2';par3';par4'];
 
%THE INTERSECTION POINT
J=[1;1;1;1];
THEINTERSECTIONPOINT =inv(Q4'*Q4)*Q4'*J
T=THEINTERSECTIONPOINT;
plot3(T(1,1),T(2,1),T(3,1),'y+')