The 23rd line from the bottom is where i attempted to write a loop. Take a look please and see if you can tell where I am going wrong. Maybe it will help show what I am trying to do here also.

Thanks,

%****************FIRST SCAN****************
%load the spreadsheet
points_1=load('C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_1_00.asc');
x_1=points_1(:,1);
y_1=points_1(:,2);
z_1=points_1(:,3);
%find centroid
N_1=size(points_1,1);
X_1=sum(x_1)/N_1;
Y_1=sum(y_1)/N_1;
Z_1=sum(z_1)/N_1;
centroid_1=[X_1 Y_1 Z_1];
%form basis for matrix elements
h1_1=x_1-X_1; h2_1=y_1-Y_1; h3_1=z_1-Z_1;

%form the elements
k11_1=sum(h1_1.*h1_1); k12_1=sum(h1_1.*h2_1); k13_1=sum(h1_1.*h3_1);
k21_1=sum(h1_1.*h2_1); k22_1=sum(h2_1.*h2_1); k23_1=sum(h2_1.*h3_1);
k31_1=sum(h1_1.*h3_1); k32_1=sum(h2_1.*h3_1); k33_1=sum(h3_1.*h3_1);

%form the matrix
M_1=[k11_1 k12_1 k13_1; k21_1 k22_1 k23_1; k31_1 k32_1 k33_1];

%find the direction vector and value
[vec_1,val_1]=eig(M_1);

%the first column of vec is the unit normal vector and corresponds to the
%minimum eigen value
un_vec_1=vec_1(:,1);
A_1=un_vec_1(1)
B_1=un_vec_1(2)
C_1=un_vec_1(3)

%calculate D using un_vec and centroid
D_1=-A_1.*centroid_1(1)-B_1.*centroid_1(2)-C_1.*centroid_1(3);

%formula for plane is A*x+B*y+C*z+D=0
%find the z's for the 4 corners of the plane
x_hi_1=max(x_1); x_lo_1=min(x_1);
y_hi_1=max(y_1); y_lo_1=min(y_1);
%z=(-D-A*x-B*y)/C
zc1_1=(-D_1-A_1*x_hi_1-B_1*y_hi_1)/C_1;
zc2_1=(-D_1-A_1*x_hi_1-B_1*y_lo_1)/C_1;
zc3_1=(-D_1-A_1*x_lo_1-B_1*y_lo_1)/C_1;
zc4_1=(-D_1-A_1*x_lo_1-B_1*y_hi_1)/C_1;

%the corner points are
cnr1_1=[x_hi_1 y_hi_1 zc1_1];
cnr2_1=[x_hi_1 y_lo_1 zc2_1];
cnr3_1=[x_lo_1 y_lo_1 zc3_1];
cnr4_1=[x_lo_1 y_hi_1 zc4_1];
p_cnrs_1=[cnr1_1; cnr2_1; cnr3_1; cnr4_1; cnr1_1];

%distance from points to plane
%use formula (Ax+By+Cz+D)/sqrt(A^2+B^2+C^2)
%the denominator = 1, just evaluate the numerator
plane2point_1=[x_1 y_1 z_1 ones(size(x_1,1),1)]*[A_1 B_1 C_1 D_1]'
m_1=1;
L=128
for n_1>0 : size(plane2point_1,1);

if L>.25
if abs(plane2point_1(n_1))<L
cull1_1(m_1,=[x_1(n_1) y_1(n_1) z_1(n_1)];
m_1=m_1+1;
L=L/2
end
end
end



%Save Point Cloud to File
fid_1=fopen('C:\MATLAB701\work\Point_Clouds1','w')
fprintf(fid_1,'%6.3f %6.3f %6.3f\n',cull1_1');
fclose(fid_1)

%Save A B C values from final plane equation
fid_1V=fopen('C:\MATLAB701\work\Vector_ABC_Values_Plane1','w')
fprintf(fid_1V,'%3.100f %3.100f %3.100f',[A_1 B_1 C_1]');
fclose(fid_1V)

%plot some points
plot3(cull1_1(:,1),cull1_1(:,2),cull1_1(:,3),'.b','MarkerSize',15);
hold on
plot3(p_cnrs_1(:,1), p_cnrs_1(:,2), p_cnrs_1(:,3),'m');
plot3(X_1,Y_1,Z_1,'.g','MarkerSize',15);