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)

Try something like:

RonL

__________________
Truth does not change because it is, or is not, believed by a majority of the people.

Giordano Bruno

Where would I place this code at in my program???

Where it says "%do calculations"

What the code fragment I posted does is loop 20 times around opening the input file writing an output file.

First file opened is the original file, subsequently it is the file just writen.

RonL

__________________
Truth does not change because it is, or is not, believed by a majority of the people.

Giordano Bruno

Where would i denote in the code to have the new files be cut down?? Maybe i need to be a bit clearer in what the code needs to do.

I import a file which contains about 15000 x y z coordinates for a point cloud. I then use this file and run my code to fit a plane to these points and tell me the distance of each point from teh plane. On the second run I take the points and delete any point that is over 128mm from the plane surface, fit a new plane to the new set of points, get the distances, and save the new point cloud file. I then run the program again with the new point cloud file and fit a newer plane to the new points, delete any points over 64mm from the surface of the plane, get distances of newer points, and save the file. I repeat this process untill i get down to deleting points .25mm away. Following is the points in my code where i delete points. Highlighted in red is the spot where the number is changed each time to tell where the cut off point for the points should be:

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

If you look at the snippet you will see that each time around the loop a different file is read and written.

The first file read from is your initial file, subsequently the file read is that written on the previous loop trip.

The names of the files written are dynamically constructed to have a numeric at the end that increases by 1 each time around the loop.

So the files written are:

'C:\MATLAB701\work\Point_Clouds11'
'C:\MATLAB701\work\Point_Clouds12'
'C:\MATLAB701\work\Point_Clouds13'

:
:

'C:\MATLAB701\work\Point_Clouds30'

RonL

__________________
Truth does not change because it is, or is not, believed by a majority of the people.

Giordano Bruno

Sorry for the confusion and being unclear on what i'm after. I not very code literate with my lingo and skills. I know in my mind what I want the code to do. Let me try to make this clearer what i'm after.

How is the points cloud being cut down each time though, this is the part where i need the loop at. I start with 128 then decrease to 64..32..16..8...4..2..1.. .5... .25. My thoughts are that you would set an integer at this part in the code (where the red 64 is at on previous reply). So set this as L. So what i want on the first run I need the value set at 128, and on the second run set L equal to 64. Basically I need a loop to decrease L by 1.5 each time untill L equals .25. To add to this I would like easily be able to change the fraction that i decrease the point cloud by for future use.

Thanks, for the help so far.

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);

You initialise a variable to 128 outside the loop, after the processing in the loop (that's in the loop at the bottom) you halve the variable, then on the next trip around the loop you have a value of 64, then 32, ...

You will find this all easier to do if you move all the actual processing inside the loop into a seperate function, so you can differentiate between the loop logic and the processing itself.

Also do not use variable names which depent on the loop index (that is don't use x_1 for the value of x first time around x_2 second time around etc. If these are not to be persistent just use x. This means you can actualy use a loop, which you can't do otherwise.

also consider designing with the aide of psuedo-code:

RonL

__________________
Truth does not change because it is, or is not, believed by a majority of the people.

Giordano Bruno

Thanks, looks like exactly what i was after at a quick glance. I wont have time to try it out till Monday morning though. Thanks for putting up with me so far. I might have a few more questions arise monday morning/late monday morning.

I hav put the code in my program, portion i added is below, but i the following error:


??? Error: File: Loop_NoP_V2_Theta_Plane1_Plane2.m Line: 19 Column: 3
Illegal use of reserved keyword "else".

Code I put in:

L=128

for idx =1 : 10
if idx==1
%load the spreadsheet
flin='C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_1_00.asc';
flout=[f100,num2str(101);
else
flin=flout
flout=[f100,num2str(100+idx)];
end

points=load(flin);


Any new helpful hints?

Thanks

Nevermind on my last post, i actaully was able to figure one out on my own, I am on the right track now. Thanks so much for you help!!

The code runs fine but only for the intial and second time through. Do you have any idea why this may be?? Code i inserted to run loop follows:

L=128
f100='C:\MATLAB701\work\Points_CloudsL1';
for idx =1 : 10
if idx==1
%load the spreadsheet
flin='C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_1_00.asc';
flout=[f100,num2str(101)];
else
flin=flout
flout=[f100,num2str(100+idx)];
end

points=load(flin);

Matlab can be flakey about some uses of the contol variable inside the loop.

Try:

L=128

for idx =1 : 10
tst=idx;
if tst==1
%load the spreadsheet
flin='C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_1_00.asc';
flout=[f100,num2str(101);
else
flin=flout
flout=[f100,num2str(100+idx)];
end

points=load(flin);

RonL

__________________
Truth does not change because it is, or is not, believed by a majority of the people.

Giordano Bruno

I have no errors when runnign my code, but it seems as if the code gives the correct values for the first two runs. I know this because I have compared outputs with the previous code I had with no "loop" in it. My code is below, please take a look and see if you can see where I need a change.

Thanks a lot.



%****************FIRST SCAN****************
L=128
f100='C:\MATLAB701\work\Points_CloudsL1';
for idx =1 : 10
if idx==1
%load the spreadsheet
flin='C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_1_00.asc';
flout=[f100,num2str(101)];
else
flin=flout
flout=[f100,num2str(100+idx)];
end

points=load(flin);


**CALCULATIONS MADE HERE NOT CURRENTLY SHOWN**

for n=1 : size(plane2point,1);
if abs(plane2point(n))<L
cull1(m,=[x(n) y(n) z(n)];
m=m+1;
end


end


L=L/2




%Save Point Cloud to File
fid_out=fopen(flout,'w')
fprintf(fid_out,'%6.3f %6.3f %6.3f\n',cull1');
fclose(fid_out)

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

end

%plot some points
plot3(cull1(:,1),cull1(:,2),cull1(:,3),'.b','MarkerSize',15);
hold on
plot3(p_cnrs(:,1), p_cnrs(:,2), p_cnrs(:,3),'m');
plot3(X,Y,Z,'.g','MarkerSize',15);

%****************SECOND SCAN****************

LI=128
f200='C:\MATLAB701\work\Points_CloudsL2';
for idxI =1 : 10
if idxI==1
%load the spreadsheet
flinI='C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_2_.02.asc';
floutI=[f200,num2str(201)];
else
flinI=floutI
floutI=[f200,num2str(200+idxI)];
end

pointsI=load(flinI);


**CALCULATIONS MADE HERE NOT CURRENTLY SHOWN**

for nI=1 : size(plane2pointI,1);
if abs(plane2pointI(nI))<LI
cull1I(mI,=[xI(nI) yI(nI) zI(nI)];
mI=mI+1;
end


end


LI=LI/2




%Save Point Cloud to File
fid_outI=fopen(floutI,'w')
fprintf(fid_outI,'%6.3f %6.3f %6.3f\n',cull1I');
fclose(fid_outI)

%Save A B C values from final plane equation
fidI=fopen('C:\MATLAB701\work\Vector_ABC_Values_Plane2L','w')
fprintf(fidI,'%3.100f %3.100f %3.100f',[AI BI CI]');
fclose(fidI)

end

%plot some points
plot3(cull1I(:,1),cull1I(:,2),cull1I(:,3),'.r','MarkerSize',15);
hold on
plot3(p_cnrsI(:,1), p_cnrsI(:,2), p_cnrsI(:,3),'m');
plot3(XI,YI,ZI,'.g','MarkerSize',15);


%**********Solving for angle between two planes***********


%Load A,B,C Values for both Planes
val_1=load('C:\MATLAB701\work\Vector_ABC_Values_Plane1L');
A1=val_1(:,1);
B1=val_1(:,2);
C1=val_1(:,3);

val_2=load('C:\MATLAB701\work\Vector_ABC_Values_Plane2L');
A2=val_2(:,1);
B2=val_2(:,2);
C2=val_2(:,3);

%Determine Magnitude of Vectors
magV1=sqrt(A1^2+B1^2+C1^2);
magV2=sqrt(A2^2+B2^2+C2^2);

M=magV1*magV2;

%Create Vectors
V1=[A1 B1 C1];
V2=[A2 B2 C2];

%Solve for angle Theta
Theta=acos((V1*V2')/M)*(180/pi)