- R2016b or later:

- R2016a or earlier (or even later, if you prefer the code)

21 views (last 30 days)

Show older comments

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?

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.

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);

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

Start Hunting!