Page 2 of 2 FirstFirst 12
Results 16 to 17 of 17

Math Help - Matlab Code Writting Help Needed

  1. #16
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by hunter View Post
    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
    Follow Math Help Forum on Facebook and Google+

  2. #17
    Newbie
    Joined
    Jul 2008
    Posts
    11

    New observation may be helpful

    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)
    Last edited by hunter; July 30th 2008 at 10:07 AM.
    Follow Math Help Forum on Facebook and Google+

Page 2 of 2 FirstFirst 12

Similar Math Help Forum Discussions

  1. Matlab code
    Posted in the Math Software Forum
    Replies: 2
    Last Post: April 4th 2011, 03:18 AM
  2. Matlab help writting an equation
    Posted in the Math Software Forum
    Replies: 1
    Last Post: December 11th 2009, 03:28 PM
  3. MATLAB code
    Posted in the Calculus Forum
    Replies: 0
    Last Post: November 5th 2009, 06:44 AM
  4. matlab code helpppp :(
    Posted in the Math Software Forum
    Replies: 3
    Last Post: October 20th 2009, 01:06 AM
  5. Need MAtlab code
    Posted in the Math Software Forum
    Replies: 1
    Last Post: April 26th 2009, 10:57 PM

Search Tags


/mathhelpforum @mathhelpforum