2
$\begingroup$

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?

$\endgroup$
6
  • $\begingroup$ The four points you’ve listed aren’t coplanar. $\endgroup$
    – amd
    Jun 30, 2017 at 0:42
  • $\begingroup$ I know. This is exactly the problem $\endgroup$
    – Davide
    Jun 30, 2017 at 0:43
  • 2
    $\begingroup$ How do you define “through,” then? $\endgroup$
    – amd
    Jun 30, 2017 at 0:44
  • $\begingroup$ Enter from a side and exit from another side. Imagine that you have a pencil and you insert into this polygon. If it pass through the polygon, then it doesn't fall to the ground. $\endgroup$
    – Davide
    Jun 30, 2017 at 0:46
  • 1
    $\begingroup$ @amd's idea is nice. Extrude the polygon to a cylinder along the direction parallel to the line. Then cut that cylinder and the line by a plane orthogonal to the line to get a polygon ans a single point. You now have reduced the problem to the well-known one of checking whether a point lies within a 2-dimensional polygon, for which there are many methods. It may be beneficial to find a 3d rotation/translation that maps your line to the z-axis. Then you can just ignore all the z-coordinates, and treat it as a 2-d problem in the xy-plane. $\endgroup$ Jun 30, 2017 at 12:56

1 Answer 1

2
$\begingroup$

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
}
$\endgroup$
8
  • $\begingroup$ Nicely done. Note that in three dimensions, G-S, even in reduced form, isn’t necessary for finding the rotated basis. The usual way I’ve seen to compute $\hat e_x$ is to take $\vec n\times\vec v_1$ and normalize. The three possible cross products with the standard basis vectors can be read off from the components of $v_1$, so there’s not really any multiplication involved this way. It’s a small, but often worthwhile, optimization. $\endgroup$
    – amd
    Jun 30, 2017 at 20:43
  • $\begingroup$ @amd: Above, we only have $\hat{e}_z$, and no hints as to the other two basis vectors, because the orientation (around final $z$ axis) is irrelevant. I use $\vec{n}$ to pick between $(1,0,0)\times\hat{e}_z$, $(0,1,0)\times\hat{e}_z$, and $(0,0,1)\times\hat{e}_z$, depending on which of $\lvert(1,0,0)\cdot\hat{e}_z\rvert$, $\lvert(0,1,0)\cdot\hat{e}_z\rvert$, $\lvert(0,0,1)\cdot\hat{e}_z\rvert$ is smallest, respectively; by simply checking which coordinate in $\hat{e}_z$ (or equivalently in $\vec{v}_1$) is the smallest in magnitude. Is that it, or is further (algorithmic) optimization possible? $\endgroup$ Jun 30, 2017 at 20:57
  • $\begingroup$ @amd: After your comment, I realized I had used $\vec{n}$ in two separate senses; as the initial axis-aligned basis vector, and as the unnormalized vector perpendicular to $\hat{e}_z$. I split those into $\vec{b}$ and $\vec{n}$, respectively; only the paragraph mentioning G-S has changed. I do believe this is the optimization you mentioned? (Also, a big thank you for looking at the math; my math-fu is pretty weak, and I do make a lot of mistakes, so even cursory checking to catch the most egregious errors is hugely appreciated!) $\endgroup$ Jun 30, 2017 at 21:04
  • $\begingroup$ What I meant what that once you’ve selected $\vec b$, you can take $\hat e_x=(\vec b\times\vec v_1)/\|\vec b\times\vec v_1\|$. Since $\vec b$ is one of the standard basis vectors, computing this cross product amounts to zeroing out the corresponding component of $\vec v_1$ and swapping the remaining two (with a change of sign). E.g., $(1,0,0)\times(x,y,z)=(0,-z,y)$. $\endgroup$
    – amd
    Jun 30, 2017 at 21:12
  • $\begingroup$ @amd: I'm starting to think I am slow.. :) But right, I now understand. Above, I calculate $\vec{n} = \vec{b} - \hat{e}_z (\hat{e}_z\cdot\vec{b})$. What you suggest yields a different $\vec{n}$, $\vec{b} \times \hat{e}_z$, which is also perpendicular to $\hat{e}_z$. It is a sensible optimization that I shall edit into my answer; thanks once again! $\endgroup$ Jun 30, 2017 at 22:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.