For loop of distance equation returning a confusing output

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!

Answer 1

[Andrew Newell]

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?

[Michael J]

This was very helpful, thank you! I am trying to keep each point of blue within a sphere of some radius around any point on the boundary of red. I don't really need to keep duplicates, so I will try and add something to the code to get rid of them. I am now running in to the issue of the code taking a long time as the matrices boundred and blue are rather large, is there an immediately apparent way to cut down on the computing required?

[Andrew Newell]

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]

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.