Using logical indexing to get subset of matrix rows which give desired output from function

11 views (last 30 days)
I have an nx3 matrix of n 3D points. I'm trying to find out how I can quickly isolate the rows of the matrix which satisfy:
dot(P1-P,normal) == 0
which would show if point P(from a given row) was in the plane given by P1 and the normal vector. I've been trying to do this using logical indexing, like
b(b>1)
in order to avoid for looping through all of the points in my space. I can't figure out how to make the logical statement dot(P1-P,normal) == 0 (instead of something like b>1), where P is a given row of my matrix. Is there a way to do this? If not, is there any way to accomplish my goal without using for loops?

Accepted Answer

Walter Roberson
Walter Roberson on 7 May 2018
Presuming that normal is a row vector (not a column vector) then for vector P1 and array P whose rows are to be dotted independently:
  • R2016b or later:
P_in_the_plane = P( conj(P1-P) * normal.' == 0, : );
  • R2016a or earlier (or even later, if you prefer the code)
P_in_the_plane = P( conj( bsxfun(@minus, P1, P) ) * normal.' == 0, : );
If your values are guaranteed real-valued then you can omit the conj()
Reminder, though, that due to round-off error, a point that is in the plane might not come out as exactly 0, so you should be considering testing the abs() against some tolerance.

More Answers (1)

Wick
Wick on 7 May 2018
Edited: Wick on 7 May 2018
I'm using random numbers just to have something to work with. The trouble is 'dot' doesn't like when the matrices aren't the same size. So we'll just use 'repmat' to make 'normal' the same size as 'P.'
You can try A == 0 for the logical indexing, but you're likely to find that machine precision will make many that "should" be 0 equal to something small, like 10*eps or smaller. I've set two logical indexes. One A == 0, the otehr abs(A) < 10*eps. You can change the multiplier in front of eps until you get what you're looking for.
P = rand(50,3);
P1 = rand(1,3);
normal = rand(1,3);
A = dot(P1-P,repmat(normal,size(P,1),1),2);
index1 = find(A == 0);
index2 = find(abs(A) < 10*eps);

Categories

Find more on Matrices and Arrays 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!