For loop of distance equation returning a confusing output

1 view (last 30 days)
Hello,
I have two matrices of coordinate points for two surfaces, one inside the other, and I am trying to create a third surface by cutting the outer surface down to a certain size (get rid of excess noise). To do this I am using two for loops, for each matrix and an if statement with the Euclidean distance equation to see if the points are within a predetermined distance. As shown below:
A = [];
for pointr = boundred
for pointb = blue
if sqrt((pointb(1)-pointr(1))^2 + (pointb(2)-pointr(2))^2 + (pointb(3)-pointr(3))^2) < 100
A = pointb;
end
end
end
The output of A is just the Z coordinates of surface blue, which confuses me. To my understanding the
for pointb = blue
iterates the loop over all of the rows, so pointb should be a row of the matrix blue. I don't understand why it only outputs pointb(3), the Z coordinate. If someone could explain why my code outputs only a part of what I am hoping to get, I would greatly appreciate it!

Accepted Answer

Andrew Newell
Andrew Newell on 24 Feb 2015
Edited: Andrew Newell on 24 Feb 2015
When something is going wrong, it's a good idea to look at the intermediate results. You could use the debugger, or you could just make some things print out as you go. As a crude check of what's going on, try the following code:
boundred = 0.3*[1 1 1; 2 2 2].'
blue = 0.2*[1 1 1; 2 2 2].'
A = [];
for pointr = boundred
pointr
for pointb = blue
pointb
if sqrt((pointb(1)-pointr(1))^2 + (pointb(2)-pointr(2))^2 + (pointb(3)-pointr(3))^2) < 100
A = pointb
end
end
end
A few things should be evident:
  1. The for statement is taking one column, not row, at a time.
  2. The value of A keeps changing.
If you want to collect the values of A, you'll need to change the assignment line to
A = [A pointb];
But if you do that, you'll find that there are repeated values of pointb in A. In essence, you're keeping each point of blue if it is within a sphere of radius 10 around any point in boundred. Is that what you want?
  3 Comments
Andrew Newell
Andrew Newell on 25 Feb 2015
The biggest slowdown in many MATLAB programs is due to commands like
A = [A pointb];
because it keeps changing the size of A (see Programming Patterns: Maximizing Code Performance by Optimizing Memory Access). You should pre-allocate vectors so MATLAB doesn't have to keep allocating new memory. But how can you do that if you don't know how large A is going to be?
By looking at the problem differently. In your code, you're deciding whether to keep each point in blue. So you could create an index idx that has an entry for each point, and then operate on that:
idx = true(1,size(blue,2));
for pointb = blue
pb = repmat(pointb,1,length(boundred));
idx = idx & sum((pb-boundred).^2) < r^2;
end
A = blue(idx);
This uses logical indexing (each component is true or false) because that is faster than numerical indexing. For 1000 points, the above code is about 100 times faster than yours.
Michael J
Michael J on 27 Feb 2015
Edited: Michael J on 27 Feb 2015
I appreciate all of your help! Though, unfortunately, I am faced with the error "Error using repmat Out of memory. Type HELP MEMORY for your options." when I use your code on my actual data. I am not sure if this error would have risen regardless of my code or not. I may have to find a way of reducing the data itself without too much of a loss of information.

Sign in to comment.

More Answers (0)

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!