I am writting a code to do the following..

I inport x y z coordinate points into the program and run my code to set a plane to the points and tell me the distance from to plain to each point. Once I know the distance I delet all points that are a certain distance away (in my case 128mm). The new set of data points is then saved to a file and imported during the second run of the program and in the end cut down even further.

As you can tell I am running a loop here, however I do not know how to write an actual loop in code to make the program delet points 128, 64, 32, 16mm...ect away. I am running thsi up to 20 times so a loop would be very useful. The first two runs of my code follow:

Like I said, basically I need help writting a loop to run instead of repeating all the same code for the "SECOND RUN". At the end of each "RUN" there is an if statement where the abs(plane2point_12(n_12))<64(hylighted to red to find easier) is the only change. I start with 128 decrease is by 1/2 each run to delet points.

I would also like to be able to change the variation at later times so if I is possable it would be great if the decrease of 1/2 could easily be changed to

1/4 or another fraction.

%****************FIRST SCAN, FIRST RUN****************

%load the spreadsheet

points_11=load('C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_1_00.asc');

x_11=points_11(:,1);

y_11=points_11(:,2);

z_11=points_11(:,3);

%find centroid

N_11=size(points_11,1);

X_11=sum(x_11)/N_11;

Y_11=sum(y_11)/N_11;

Z_11=sum(z_11)/N_11;

centroid_11=[X_11 Y_11 Z_11];

%form basis for matrix elements

h1_11=x_11-X_11; h2_11=y_11-Y_11; h3_11=z_11-Z_11;

%form the elements

k11_11=sum(h1_11.*h1_11); k12_11=sum(h1_11.*h2_11); k13_11=sum(h1_11.*h3_11);

k21_11=sum(h1_11.*h2_11); k22_11=sum(h2_11.*h2_11); k23_11=sum(h2_11.*h3_11);

k31_11=sum(h1_11.*h3_11); k32_11=sum(h2_11.*h3_11); k33_11=sum(h3_11.*h3_11);

%form the matrix

M_11=[k11_11 k12_11 k13_11; k21_11 k22_11 k23_11; k31_11 k32_11 k33_11];

%find the direction vector and value

[vec_11,val_11]=eig(M_11);

%the first column of vec is the unit normal vector and corresponds to the

%minimum eigen value

un_vec_11=vec_11(:,1);

A_11=un_vec_11(1)

B_11=un_vec_11(2)

C_11=un_vec_11(3)

%calculate D using un_vec and centroid

D_11=-A_11.*centroid_11(1)-B_11.*centroid_11(2)-C_11.*centroid_11(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_11=max(x_11); x_lo_11=min(x_11);

y_hi_11=max(y_11); y_lo_11=min(y_11);

%z=(-D-A*x-B*y)/C

zc1_11=(-D_11-A_11*x_hi_11-B_11*y_hi_11)/C_11;

zc2_11=(-D_11-A_11*x_hi_11-B_11*y_lo_11)/C_11;

zc3_11=(-D_11-A_11*x_lo_11-B_11*y_lo_11)/C_11;

zc4_11=(-D_11-A_11*x_lo_11-B_11*y_hi_11)/C_11;

%the corner points are

cnr1_11=[x_hi_11 y_hi_11 zc1_11];

cnr2_11=[x_hi_11 y_lo_11 zc2_11];

cnr3_11=[x_lo_11 y_lo_11 zc3_11];

cnr4_11=[x_lo_11 y_hi_11 zc4_11];

p_cnrs_11=[cnr1_11; cnr2_11; cnr3_11; cnr4_11; cnr1_11];

%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_11=[x_11 y_11 z_11 ones(size(x_11,1),1)]*[A_11 B_11 C_11 D_11]'

m_11=1;

for n_11=1 : size(plane2point_11,1);

if abs(plane2point_11(n_11))<128

cull1_11(m_11,=[x_11(n_11) y_11(n_11) z_11(n_11)];

m_11=m_11+1;

end

end

fid_11=fopen('C:\MATLAB701\work\Point_Clouds11','w')

fprintf(fid_11,'%6.3f %6.3f %6.3f\n',cull1_11');

fclose(fid_11)

%****************FIRST SCAN,SECOND RUN****************

%load the spreadsheet

points_12=load('C:\MATLAB701\work\Point_Clouds11');

x_12=points_12(:,1);

y_12=points_12(:,2);

z_12=points_12(:,3);

%find centroid

N_12=size(points_12,1);

X_12=sum(x_12)/N_12;

Y_12=sum(y_12)/N_12;

Z_12=sum(z_12)/N_12;

centroid_12=[X_12 Y_12 Z_12];

%form basis for matrix elements

h1_12=x_12-X_12; h2_12=y_12-Y_12; h3_12=z_12-Z_12;

%form the elements

k11_12=sum(h1_12.*h1_12); k12_12=sum(h1_12.*h2_12); k13_12=sum(h1_12.*h3_12);

k21_12=sum(h1_12.*h2_12); k22_12=sum(h2_12.*h2_12); k23_12=sum(h2_12.*h3_12);

k31_12=sum(h1_12.*h3_12); k32_12=sum(h2_12.*h3_12); k33_12=sum(h3_12.*h3_12);

%form the matrix

M_12=[k11_12 k12_12 k13_12; k21_12 k22_12 k23_12; k31_12 k32_12 k33_12];

%find the direction vector and value

[vec_12,val_12]=eig(M_12);

%the first column of vec is the unit normal vector and corresponds to the

%minimum eigen value

un_vec_12=vec_12(:,1);

A_12=un_vec_12(1)

B_12=un_vec_12(2)

C_12=un_vec_12(3)

%calculate D using un_vec and centroid

D_12=-A_12.*centroid_12(1)-B_12.*centroid_12(2)-C_12.*centroid_12(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_12=max(x_12); x_lo_12=min(x_12);

y_hi_12=max(y_12); y_lo_12=min(y_12);

%z=(-D-A*x-B*y)/C

zc1_12=(-D_12-A_12*x_hi_12-B_12*y_hi_12)/C_12;

zc2_12=(-D_12-A_12*x_hi_12-B_12*y_lo_12)/C_12;

zc3_12=(-D_12-A_12*x_lo_12-B_12*y_lo_12)/C_12;

zc4_12=(-D_12-A_12*x_lo_12-B_12*y_hi_12)/C_12;

%the corner points are

cnr1_12=[x_hi_12 y_hi_12 zc1_12];

cnr2_12=[x_hi_12 y_lo_12 zc2_12];

cnr3_12=[x_lo_12 y_lo_12 zc3_12];

cnr4_12=[x_lo_12 y_hi_12 zc4_12];

p_cnrs_12=[cnr1_12; cnr2_12; cnr3_12; cnr4_12; cnr1_12];

%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_12=[x_12 y_12 z_12 ones(size(x_12,1),1)]*[A_12 B_12 C_12 D_12]'

m_12=1;

for n_12=1 : size(plane2point_12,1);

if abs(plane2point_12(n_12))<64

cull1_12(m_12,=[x_12(n_12) y_12(n_12) z_12(n_12)];

m_12=m_12+1;

end

end

fid_12=fopen('C:\MATLAB701\work\Point_Clouds12','w')

fprintf(fid_12,'%6.3f %6.3f %6.3f\n',cull1_12');

fclose(fid_12)