Thread: Matlab Code Writting Help Needed

Originally Posted by hunter

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(101+idx)];
end

points=load(flin);

See above in red, but that should not have stopped it working!

RonL

I have noticed that the output files are being decreased any. They should be since I am deleting points. No points get deleted untill after the second scan which is why the "looped" version and my "old" version outputs numbers match up. Any input on how to fix this?? If needed i can send you my "old" none looped version to take a look at also. The current code full with the loop that I am running follows:

%****************FIRST SCAN****************
L=128
f1000='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=[f1000,num2str(1001)];
else
flin=flout
flout=[f1000,num2str(1001+idx)];
end

points=load(flin);


x=points(:,1);
y=points(:,2);
z=points(:,3);
%find centroid
N=size(points,1);
X=sum(x)/N;
Y=sum(y)/N;
Z=sum(z)/N;
centroid=[X Y Z];
%form basis for matrix elements
h1=x-X; h2=y-Y; h3=z-Z;

%form the elements
k11=sum(h1.*h1); k12=sum(h1.*h2); k13=sum(h1.*h3);
k21=sum(h1.*h2); k22=sum(h2.*h2); k23=sum(h2.*h3);
k31=sum(h1.*h3); k32=sum(h2.*h3); k33=sum(h3.*h3);

%form the matrix
M=[k11 k12 k13; k21 k22 k23; k31 k32 k33];

%find the direction vector and value
[vec,val]=eig(M);

%the first column of vec is the unit normal vector and corresponds to the
%minimum eigen value
un_vec=vec(:,1);
A=un_vec(1)
B=un_vec(2)
C=un_vec(3)

%calculate D using un_vec and centroid
D=-A.*centroid(1)-B.*centroid(2)-C.*centroid(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=max(x); x_lo=min(x);
y_hi=max(y); y_lo=min(y);
%z=(-D-A*x-B*y)/C
zc1=(-D-A*x_hi-B*y_hi)/C;
zc2=(-D-A*x_hi-B*y_lo)/C;
zc3=(-D-A*x_lo-B*y_lo)/C;
zc4=(-D-A*x_lo-B*y_hi)/C;

%the corner points are
cnr1=[x_hi y_hi zc1];
cnr2=[x_hi y_lo zc2];
cnr3=[x_lo y_lo zc3];
cnr4=[x_lo y_hi zc4];
p_cnrs=[cnr1; cnr2; cnr3; cnr4; cnr1];

%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=[x y z ones(size(x,1),1)]*[A B C D]';
m=1;
for n=1 : size(plane2point,1);
if abs(plane2point(n))<L
cull1(m,=[x(n) y(n) z(n)];
m=m+1;


end

L=L/2

end



%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
f1000='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_1_00.asc';
floutI=[f1000,num2str(1001)];
else
flinI=floutI
floutI=[f1000,num2str(1001+idx)];
end

pointsI=load(flinI);

xI=pointsI(:,1);
yI=pointsI(:,2);
zI=pointsI(:,3);
%find centroid
NI=size(pointsI,1);
XI=sum(xI)/NI;
YI=sum(yI)/NI;
ZI=sum(zI)/NI;
centroidI=[XI YI ZI];
%form basis for matrix elements
h1I=xI-XI; h2I=yI-YI; h3I=zI-ZI;

%form the elements
k11I=sum(h1I.*h1I); k12I=sum(h1I.*h2I); k13I=sum(h1I.*h3I);
k21I=sum(h1I.*h2I); k22I=sum(h2I.*h2I); k23I=sum(h2I.*h3I);
k31I=sum(h1I.*h3I); k32I=sum(h2I.*h3I); k33I=sum(h3I.*h3I);

%form the matrix
MI=[k11I k12I k13I; k21I k22I k23I; k31I k32I k33I];

%find the direction vector and value
[vecI,valI]=eig(MI);

%the first column of vec is the unit normal vector and corresponds to the
%minimum eigen value
un_vecI=vecI(:,1);
AI=un_vecI(1)
BI=un_vecI(2)
CI=un_vecI(3)

%calculate D using un_vec and centroid
DI=-AI.*centroidI(1)-BI.*centroidI(2)-CI.*centroidI(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_hiI=max(xI); x_loI=min(xI);
y_hiI=max(yI); y_loI=min(yI);
%z=(-D-A*x-B*y)/C
zc1I=(-DI-AI*x_hiI-BI*y_hiI)/CI;
zc2I=(-DI-AI*x_hiI-BI*y_loI)/CI;
zc3I=(-DI-AI*x_loI-BI*y_loI)/CI;
zc4I=(-DI-AI*x_loI-BI*y_hiI)/CI;

%the corner points are
cnr1I=[x_hiI y_hiI zc1I];
cnr2I=[x_hiI y_loI zc2I];
cnr3I=[x_loI y_loI zc3I];
cnr4I=[x_loI y_hiI zc4I];
p_cnrsI=[cnr1I; cnr2I; cnr3I; cnr4I; cnr1I];

%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
plane2pointI=[xI yI zI ones(size(xI,1),1)]*[AI BI CI DI]';
mI=1;
for nI=1 : size(plane2pointI,1);
if abs(plane2pointI(nI))<LI
cull1I(mI,=[xI(nI) yI(nI) zI(nI)];
mI=mI+1;


end

LI=LI/2

end


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