Intersection of a 2D polygon in 3D and a straight line

I have a polygon determined by 5 points in 3D. I want to see if the line, say X-axis interescts with this polygon. Also, I am interested in the intersection point.

For instance, the polygon is determined by the following points that are coplanar:

V1 = [11.2908113470153 -8.47779952423983 8.18414177879808]

V2 = [12.3366023247650 -8.58127758366331 8.28004384038752]

V3 = [1.2757044711089 -8.47363504164875 8.18612938262431]

V4 = [13.3363266184718 -7.53719741280684 9.81580879134791]

V5 = [15.0742258082624 -7.67475550912392 10.0186434441811].

The end points of the line segment along the X-axis are [-30 0 0] and [30 0 0].

Answer 1

[Bjorn Gustavsson]

OK, you'll get a bit of pseudo-code/algebra, that solves most of the problem. Lets write the equation for the plane as:

that is all points in 3-D r who's scalar product with the unit-vector, normal to the plane, is l. We then write the equation for the line:

This gives us the equation for the line-plane intersection:

Which we can solve for s at the intersection:

giving us the intersection-point:

So that makes the first step.

Next step I'd take would be to make a 2-D coordinate system in the plane, any random origin-point will do (one of your polygon-points seems convenient), then a unit-vector towards one of its neighbors will make a good direction, and an arbitrary direction perpendicular to that. Now you have an origin and 2 basis-vectors for your plane and can convert your 5+1 (minus the origin-point...) to their corresponding points in this plane. That completes the second step.

Third step is to use the matlab-function inpolygon - which gives you the in/out-answer.

I leave the 2 last steps for you to do, it shold be rather straightforward, and makes for a nice geometric function, you can keep forever.

HTH

[LC]

Thanks! I think I lost you at Step 2. I got Steps 1 and 2. Can you please elaborate on Step 2?

[Bjorn Gustavsson]

Since you have a plane problem we can transform your points into a coordinate system, lets call it [u,v,w], where all 6 points of yours will have w-coordinate 0. That transformation will consist of one translation and a multiplication with a rotation-matrix. You can create the rotation-matrix from , and your other points, make the first in-plane unit-vector from the first 2 points:

and the second in-plane unit-vector you get by cross-multiplying these two:

This will give you the rotation-matrix:

(the latex-notation fell appart here, the unit-vectors should make up rows in M) Then you get your 2-D coordinates by:

all should have their third element equal to zero, so these you could then feed to inpolygon.

Another way of doing this is to look at the cross-product between the vectors and - this gives you a simple mechanical solution that works under some conditions - if you look at the simple case where you have all points already in the x-y plane you can figure out what to look for and what the condition is.

HTH

[LC]

Ok. So I tried your suggestion. I get the error at the last line of the code. I am unsure where things went wrong. I am showing the small code that I wrote here.

clc
clear
L = 30;
r1 = - [L 0 0]; % First endpoint of the line segment
r0 = [L 0 0]; % Second endpoint of the line segment
v = [-1,-1,1;-1,-1,-1;-1,1,-1;-1,1,1]; % Let us take a facet of cuboid for simplicity
vert1 = v(1,:); % Pick the first vertex to calculate the normal vector to the facet

pol_ang =  atan2(v(:,1),v(:,2));
vert_ang = [v pol_ang];
vert_sort = sortrows(vert_ang,4); % Sorting vertices of facet anticlockwise
vert_sort(:,4) = [];

%facet normal vector
norm_vec = cross(vert_sort(2,:) - vert_sort(1,:), vert_sort(length(vert_sort),:)- vert_sort(1,:));
e_n = norm_vec/norm(norm_vec);
ll = dot(v(1,:),e_n);
e_l = (r1-r0)/norm(r1-r0);
r_l = r0 + (ll - dot(r0,e_n))/(dot(e_l, e_n))*e_l;
e_1 = (v(2,:) - v(1,:))/(norm(v(2,:) - v(1,:)));
e_2 = cross(e_n,e_1);
M = [e_1;e_2;e_n];
v_new= dot(M,(v-r0));

[Bjorn Gustavsson]

OK, you're almost there.

1, I goofed with the notation. I happened to use twice. In the first answer it denotes the first point on your line. In the second I unfortunately and reclessly used it to denote one of the vertex-points in your plane.

That means you should change your last line to:

v_new= M*(v-v(1,:))';

This means we set the new origin to the first vertex. Here you should also note that the multiplication between your vertices and the rotation-matrix is a vector-matrix multiplication (or 4 of them). (This is a mathematical notation overload, that actually doesn't cause confusion - once you get used to it...) To do the vectors-matrix multiplication in one swoop you have to have the coordinates in the column-direction (this causes me frustration, because my prefered orientation of vectors like this is the one you use, there will be additional transpose-operations to pay for this choise). So that gives you the new 3-D coordinates of the vertices. You have to do the same transformation for your point-projection:

r_lnew= M*(r_l-v(1,:))';

This should set you up nicely for using inpolygon:

inpolygon(r_lnew(1),r_lnew(2),v_new(1,:),v_new(2,:))

Another advice for geometry-problems like this is to plot things along the way - I find that this helps me "anchor my thinking" so I don't spin out of order...

HTH