Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It only takes a minute to sign up.
Sign up to join this community
Is there a way to determine if a straight line passes through a polygon?
For example we have a polygon with: A=(1, 1, 1), B=(1.5, 1, 1), C=(1.7, 1.7, 1.2), D=(0.9, 1.4, 0.9)
The line (1.2, 1.1, 2)+t(0.1, 0.2, -1) pass through the polygon but the line (2, 2, 3)+t(-0.1, 1, -1.1) not.
Is there a formula?
I am shamelessly stealing the ideas expressed in comments to the question by amd and Jaap Scherphuis: we rotate and translate the coordinate system so that the line is the $z$ axis, and then find if the polygon projected to the $xy$ plane contains the origin.
My intention is only to show how a programmer might apply the suggestions to solve this problem; all kudos belongs to amd and Jaap Scherphuis.
First, let's call the vertices of the polygon $$\vec{p}_n = ( x_n , y_n , z_n ), \qquad n = 0 .. N-1$$ and define the line $\vec{v}(t)$, $t \in \mathbb{R}$ as $$\vec{v}(t) = \vec{v}_0 + t \vec{v}_1$$
In OP's case, $\vec{v}_0 = ( x_{v0} , y_{v0} , z_{v0} ) = ( 1.2 , 1.1 , 2 )$, $\vec{v}_1 = ( x_{v1} , y_{v1} , z_{v1} ) = ( 0.1 , 0.2 , -1 )$, $N = 4$, $\vec{p}_0 = ( 1, 1, 1 )$, $\vec{p}_1 = ( 1.5 , 1 , 1 )$, $\vec{p}_2 = ( 1.7, 1.7, 1.2 )$, and $\vec{p}_3 = ( 0.9, 1.4, 0.9 )$.
The translation that moves the line so that it passes through the origin is obviously $-\vec{v}_0$ (applied before the rotation). We need to find an orthonormal rotation matrix $\mathbf{M}$ that rotates the line to the $z$ axis.
There is a "trick" we can use here: orthogonal matrices' inverse is their transpose. In essence, if we find $\mathbf{M}^{-1} = \mathbf{M}^{T}$ that rotates the $z$ axis to the line, we immediately find $\mathbf{M}$. Let's define $\hat{e}_x$, $\hat{e}_y$, and $\hat{e}_z$ as the columns of $\mathbf{M}$: $$\mathbf{M} = \left [ \hat{e}_x \; \hat{e}_y \; \hat{e}_z \right ] = \left [ \begin{array}{ccc} X_x & Y_x & Z_x \\ X_y & Y_y & Z_y \\ X_z & Y_z & Z_z \end{array} \right ]$$ where$$\begin{array}{c} \hat{e}_x = ( X_x , Y_x , Z_x ) \\ \hat{e}_y = ( Y_x , Y_y , Z_y ) \\ \hat{e}_z = ( Z_x , Z_y , Z_z ) \end{array}$$
We already know the $z$ vector; we just need to normalize it to length 1: $$\hat{e}_z = \frac{\vec{v}_1}{\lVert \vec{v}_1 \rVert} = \frac{\vec{v}_1}{\sqrt{\vec{v}_1 \cdot \vec{v}_1}}$$i.e.$$\begin{cases} Z_x = \frac{x_{v1}}{\sqrt{x_{v1}^2 + y_{v1}^2 + z_{v1}^2}} \\ Z_y = \frac{y_{v1}}{\sqrt{x_{v1}^2 + y_{v1}^2 + z_{v1}^2}} \\ Z_z = \frac{z_{v1}}{\sqrt{x_{v1}^2 + y_{v1}^2 + z_{v1}^2}} \end{cases}$$
Because we don't care about the orientation of the 2D polygon -- that is, its rotation around the $z$ axis --, any two vectors $\hat{e}_x$ and $\hat{e}_y$ that are of unit length and orthogonal (to each other and to $\hat{e}_z$), will work. In fact, we only need $\hat{e}_y$, because $\mathbf{M}$ being orthonormal means $\hat{e}_x = \hat{e}_y \times \hat{e}_z$.
One easy way to find an orthogonal vector is to start with the unit axis vectors $(1, 0, 0)$, $(0, 1, 0)$, and $(0, 0, 1)$, and pick the one that is furthest away from $\hat{e}_z$. Let's say this is $\vec{b}$. We then orthonormalize it against $\hat{e}_z$ using the Gram-Schmidt process: $\vec{n} = \vec{b} - \hat{e}_z ( \hat{e}_z \cdot \vec{b} )$, $\hat{e}_y = \vec{n} / \lVert \vec{n} \rVert$. The last orthogonal vector in the basis is then $\hat{e}_x = \hat{e}_y \times \hat{e}_z$.
However, in a comment, amd suggests using $\vec{n} = \vec{b} \times \hat{e}_z$ instead; then, $\hat{e}_x = \vec{n} / \lVert \vec{n} \rVert$ and $\hat{e}_y = \hat{e}_z \times \hat{e}_x$, giving us the same basis vectors $\hat{e}_x$, $\hat{e}_y$, and $\hat{e}_z$ as above. Because $\vec{b}$ is one of the unit axis vectors, we only need to find the smallest component of $\hat{e}_z$ (or, equivalently, $\vec{v}_1$) in magnitude, zero the corresponding component, swap the two other components, and negate one of them, to get $\vec{n}$:
If $\lvert Z_x \rvert \le \lvert Z_y \rvert$ and $\lvert Z_x \rvert \le \lvert Z_z \rvert$, then
$\vec{n} = ( 0 ,\, -Z_z ,\, Z_y )$
otherwise,
If $\lvert Z_y \rvert \le \lvert Z_x \rvert$ and $\lvert Z_y \rvert \le \lvert Z_z \rvert$, then
$\vec{n} = ( Z_z ,\, 0 ,\, -Z_y )$
otherwise,
$\vec{n} = ( -Z_y ,\, Z_x ,\, 0 )$
We normalize $\vec{n}$ to get the first column vector $\hat{e}_x$: $$\hat{e}_x = \frac{\vec{n}}{\lVert \vec{n} \rVert} = \frac{\vec{n}}{\sqrt{\vec{n} \cdot \vec{n}}}$$ and the second column vector of $\mathbf{M}$, $\hat{e}_y$, is $$\hat{e}_y = \hat{e}_z \times \hat{e}_x = \left ( X_z Z_y - X_y Z_z ,\; X_x Z_z - X_z Z_x ,\; X_y Z_x - X_x Z_y \right )$$ (Although mathematically $\hat{e}_y$ will have unit length, you may wish to normalize it to unit length by dividing each component by $\sqrt{Y_x^2 + Y_y^2 + Y_z^2}$ in case of numerical errors, like rounding errors. Normally, this is not needed.)
Now that we know the rotation matrix $\mathbf{M}$ and the translation before rotation $-\vec{v}_0$, we can transform and project the polygon vertices to the $z = 0$ plane. Essentially, $${\vec{p}'}_n = \mathbf{M} \left ( \vec{p}_n - \vec{v}_0 \right )$$ Component-wise, this boils down to $$\begin{array}{c} {x'}_n = X_x ( x_n - x_{v0} ) + X_y ( y_n - y_{v0} ) + X_z ( z_n - z_{v0} ) \\ {y'}_n = Y_x ( x_n - x_{v0} ) + Y_y ( y_n - y_{v0} ) + Y_z ( Z_n - z_{v0} ) \end{array}$$ We do not need the ${z'}_n$ coordinates, because projecting to the $z = 0$ plane sets all $z$ components to zero.
At this point, the problem has been simplified to checking if the origin $(0, 0)$ is inside the 2D polygon defined by vertices $({x'}_n , {y'}_n )$: the well-known point in polygon test. There are two basic approaches: the ray casting test, and the winding number test.
In this particular case, I warmly recommend the winding number test, using quadrants or octants (or octants plus axes), i.e. 4, 8, or 12 directions only, using only coordinate comparisons, no trigonometric functions at all.
The reason is that this way, you don't need a separate array for the ${x'}_n$ and ${y'}_n$ coordinates: you simply check which direction (quadrant/octant) the point $({x'}_n , {y'}_n)$ falls, and update the winding accordingly. (You do need to remember the previous point, so that you can solve the cases where the polygon edge skips to the opposite side, to find out which side (clockwise or counterclockwise) it went.)
Here is an example awk script, that takes six parameters (the components of $\vec{v}_0$ and $\vec{v}_1$), applies the translation and rotation to the 3D points specified in standard input, one point per line, and projects the results to the $z=0$ plane, printing the 2D points to standard output:
#!/usr/bin/awk -f
function abs(value) { if (value < 0) return -value; return value }
BEGIN {
if (":" epsilon == ":")
epsilon = 0.000005;
if (length(ARGV) != 7) {
printf "\n" > "/dev/stderr"
printf "Usage: %s [ -h | --help ]\n", ARGV0 > "/dev/stderr"
printf " %s X0 Y0 Z0 X1 Y1 Z1 < FILE\n", ARGV0 > "/dev/stderr"
printf "\n" > "/dev/stderr"
printf "This script reads 3D points from FILE,\n" > "/dev/stderr"
printf "translating them so that X0,Y0,Z0 moves to origin,\n" > "/dev/stderr"
printf "rotates them so that X1,Y1,Z1 is on the positive z axis,\n" > "/dev/stderr"
printf "projects the transformed points to the xy plane,\n" > "/dev/stderr"
printf "and outputs the resulting 2D coordinates.\n" > "/dev/stderr"
printf "\n" > "/dev/stderr"
exit(1)
}
ARGV0 = ARGV[0]
ARGV1 = ARGV[1]
ARGV2 = ARGV[2]
ARGV3 = ARGV[3]
ARGV4 = ARGV[4]
ARGV5 = ARGV[5]
ARGV6 = ARGV[6]
split("", ARGV)
ARGV[0] = ARGV0
X0 = ARGV1 - 0
Y0 = ARGV2 - 0
Z0 = ARGV3 - 0
Zx = ARGV4 - 0
Zy = ARGV5 - 0
Zz = ARGV6 - 0
Zn = sqrt(Zx*Zx + Zy*Zy + Zz*Zz)
if (Zn <= epsilon) {
printf "%s %s %s: Vector is too short.\n", ARGV4, ARGV5, ARGV6 > "/dev/stderr"
exit(1)
}
Zx = Zx / Zn
Zy = Zy / Zn
Zz = Zz / Zn
if (abs(Zx) <= abs(Zy) && abs(Zx) <= abs(Zz)) {
Xx = 0
Xy = -Zz
Xz = Zy
} else
if (abs(Zy) <= abs(Zx) && abs(Zy) <= abs(Zz)) {
Xx = Zz
Xy = 0
Xz = -Zx
} else {
Xx = -Zy
Xy = Zx
Xz = 0
}
Xn = sqrt(Xx*Xx + Xy*Xy + Xz*Xz)
if (Xn <= epsilon) {
printf "Oops: X basis vector became zero length.\n" > "/dev/stderr"
exit(1)
}
Xx = Xx / Xn
Xy = Xy / Xn
Xz = Xz / Xn
Yx = Xz*Zy - Xy*Zz
Yy = Xx*Zz - Xz*Zx
Yz = Xy*Zx - Xx*Zy
Yn = sqrt(Yx*Yx + Yy*Yy + Yz*Yz)
printf "T = %11.6f %11.6f %11.6f\n", -X0, -Y0, -Z0 > "/dev/stderr"
printf "eX = %11.6f %11.6f %11.6f\n", Xx, Xy, Xz > "/dev/stderr"
printf "eY = %11.6f %11.6f %11.6f\n", Yx, Yy, Yz > "/dev/stderr"
printf "eZ = %11.6f %11.6f %11.6f\n", Zx, Zy, Zz > "/dev/stderr"
printf "eX.eX = %8.6f eX.eY = %8.6f\n", Xx*Xx+Xy*Xy+Xz*Xz, Xx*Yx+Xy*Yy+Xz*Yz > "/dev/stderr"
printf "eY.eY = %8.6f eX.eZ = %8.6f\n", Yx*Yx+Yy*Yy+Yz*Yz, Xx*Zx+Xy*Zy+Xz*Zz > "/dev/stderr"
printf "eZ.eZ = %8.6f eY.eZ = %8.6f\n", Zx*Zx+Zy*Zy+Zz*Zz, Yx*Zx+Yy*Zy+Yz*Zz > "/dev/stderr"
}
NF >= 3 {
x = $1 - X0
y = $2 - Y0
z = $3 - Z0
printf "%16.6f %16.6f\n", Xx*x + Xy*y + Xz*z, Yx*x + Yy*y + Yz*z
}
