Page 1 of 2 12 LastLast
Results 1 to 15 of 17

Math Help - Matlab Code Writting Help Needed

  1. #1
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Matlab Code Writting Help Needed

    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)
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Try something like:

    Code:
    fl0='C:\MATLAB701\work\Point_Clouds';
    for idx =1 : 20
      if idx==1
        flin='C:\Documents and Settings\prhill\My Documents\Hill\Point_Clouds\20_1_00.asc';
        flout=[fl0,num2str(11);
      else
        flin=flout
        flout=[fl0,num2str(10+idx)];
      end
      
      points=load(flin);
      
      %do calcuations
      
      fid_out=fopen(flout,'w');
      fprintf(fid_out,'%6.3f %6.3f %6.3f\n',cull1');
      fclose(fid_out);
      
    end
    RonL
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Placement of code

    Where would I place this code at in my program???
    Follow Math Help Forum on Facebook and Google+

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

  5. #5
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Sounds good, but i think i need a bit more if you don't mind

    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
    Follow Math Help Forum on Facebook and Google+

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

  7. #7
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Understand that part, but where do you denote the decrease in point cloud

    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.
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Here is my failed loop attempt

    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);
    Follow Math Help Forum on Facebook and Google+

  9. #9
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by hunter View Post
    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.
    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:

    Code:
    L=128
    
    for idx =1 : 20
      if idx==1
        read from original input file
        set file name for first output file
      else
        read from last output file
        set file name for next output file
      end
      
    
      
      %do calcuations
      
      write output file
      do any other plotting, printing etc
      L=L/2
      
    end
    RonL
    Follow Math Help Forum on Facebook and Google+

  10. #10
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Looks Good

    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.
    Follow Math Help Forum on Facebook and Google+

  11. #11
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Tried Code, but error received

    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
    Follow Math Help Forum on Facebook and Google+

  12. #12
    Newbie
    Joined
    Jul 2008
    Posts
    11

    Disregard Last Post

    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!!
    Follow Math Help Forum on Facebook and Google+

  13. #13
    Newbie
    Joined
    Jul 2008
    Posts
    11

    New Question

    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);
    Follow Math Help Forum on Facebook and Google+

  14. #14
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by hunter View Post
    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

    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
    Follow Math Help Forum on Facebook and Google+

  15. #15
    Newbie
    Joined
    Jul 2008
    Posts
    11

    New problem

    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)

    Follow Math Help Forum on Facebook and Google+

Page 1 of 2 12 LastLast

Similar Math Help Forum Discussions

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

Search Tags


/mathhelpforum @mathhelpforum