Way to check if a 3D polyline fully going around an axis?

Question

I have an application where it is important to check if a 3D polyline revolve fully around an axis. Here's two example where the polylines do revolve around Oz:

And here is an example where it does not:

One constraint I have that all polylines are connected (i.e the first and the last vertex is the same). Any ideas? Thanks.

Answer 1

Another option is to rotate and translate the polyline so that the axis passes through origin, and is parallel to the $z$ axis. Then, this simplifies to the classic 2D point-in-polygon test.

I personally prefer the winding number check using "octants": $$\begin{array}{c|c} \text{Rule} & \text{If True} \\ \hline \Delta x \gt 0 & +1 \\ \Delta x \lt 0 & +2 \\ \Delta y \gt 0 & +3 \\ \Delta y \lt 0 & +6 \\ \lvert \Delta x \rvert \lt \lvert \Delta y \rvert & +9 \\ \lvert \Delta x \rvert \gt \lvert \Delta y \rvert & +18 \\ \end{array} \quad \begin{array}{c|c} \text{Sum} & \text{Octant} \\ \hline 19 & 0 \\ 4 & 1 \\ 12 & 2 \\ 5 & 3 \\ 20 & 4 \\ 8 & 5 \\ 15 & 6 \\ 7 & 7 \\ \end{array} ~ \begin{array}{c|c} \text{Sum} & \text{Octant} \\ \hline 22 & 0 \\ 13 & 1 \\ 14 & 2 \\ 23 & 3 \\ 26 & 4 \\ 17 & 5 \\ 16 & 6 \\ 25 & 7 \\ \end{array}$$ Sum $0$ indicates origin, and all other sums are impossible. Here is an illustration of the octants and the sums using the above rules: Note that odd octants include a diagonal, and even octants an axis. This is one way to ensure the below winding counting handles nontrivial cases correctly.

Let $o_i$ be the octant for vertex $i$ in the polyline with respect to the axis, for $i = 0, \dots\ n-1$, and for simplicity, $o_n = o_0$ for the repeated closing vertex. Then, the winding number $w$ for the polyline is $$w = \sum_{i=0}^{n-1} (((12 + o_{i+1} - o_{i}) \mod 8) - 4)$$ except for the cases where $o_{i+1} - o_{i} = \pm 4$, which are handled separately:

• Calculate $$\begin{aligned} \Delta x &= x_{i+1} - x_{i} \\ \Delta y &= y_{i+1} - y_{i} \\ \end{aligned}$$ noting that $x_{i}$, $y_{i}$, $x_{i+1}$, and $y_{i+1}$ do not necessarily need to have the axis at origin, as long as they are Cartesian coordinates in a plane perpendicular to the axis. Usually they'll be with the axis at origin, though.

• If $\lvert \Delta x \lvert \ge \lvert \Delta y \rvert$, then calculate $$\delta y = y_{i} x_{i+1} - x_{i} y_{i+1}$$ (The line segment intersects $x = 0$ at $y = \delta y / \Delta x$.)
  If $\delta y = 0$, the polyline passes through the axis.
  If $\delta y \gt 0$, subtract $4$ from $w$.
  If $\delta y \lt 0$, add $4$ to $w$.

• Otherwise, calculate $$\delta x = x_{i} y_{i+1} - y_{i} x_{i+1}$$ (The line segment intersects $y = 0$ at $x = \delta x / \Delta y$.)
  If $\delta x = 0$, the polyline passes through the axis.
  If $\delta x \gt 0$, add $4$ to $w$.
  If $\delta x \lt 0$, subtract $4$ from $w$.

The expression $((12 + o_{i+1} - o_{i})\mod 8) - 4$ just maps the difference to range $-4$ to $+3$, inclusive: $$\begin{array}{c|c} \text{difference} & \text{value} \\ \hline -7 & 1 \\ -6 & 2 \\ -5 & 3 \\ -4 & \text{special} \\ -3 & -3 \\ -2 & -2 \\ -1 & -1 \\ 0 & 0 \\ 1 & 1 \\ 2 & 2 \\ 3 & 3 \\ 4 & \text{special} \\ 5 & -3 \\ 6 & -2 \\ 7 & -1 \\ \end{array}$$ In Python and C and many other languages, with o0 as $o_i$ and o1 as $o_{i+1}$, you can use ((12 + o0 - o1) & 7) - 4.

For a closed polyline, $w = 8 k$, where $k \in \mathbb{Z}$.
If $w = 0$, axis is outside the polyline.
If $w = \pm 8$, axis in inside the polyline.
Otherwise, $w = \pm 16, \pm 24, \pm 32, \dots$, the polyline loops around the axis more than once. Depending on the situation, you probably wish to consider the axis "inside" the polyline in these cases. (In 2D, there is a choice: nonzero considers the axis to be inside the polyline in those cases; odd-even considers the axis to be inside the polyline if $w/8$ is odd, and outside if even. But this does not really apply to the 3D polyline.)

I prefer this approach over the intersection counting approach, because this tends to be faster for longer polylines (more vertices) and easier to implement in a numerically robust way, as there is only one "special case" (when a line segment passes through four octants), and that too has a simple solution with just additions, subtractions, multiplications, and comparisons.

---

Here is an example Python 3 implementation:

from math import sqrt

def octant(dx, dy):
    """Origin is None, and the octants are numbered 0 to 7,
       starting from the positive X axis towards the positive Y axis.
       Even octants include an axis, odd octants a diagonal."""
    return (None,None,None,None,1,3,None,7,5,
            None,None,None,   2,1,2,   6,6,5,
            None,   0,   4,None,0,3,None,7,4)[ 1*(dx > 0)            +
                                               2*(dx < 0)            +
                                               3*(dy > 0)            +
                                               6*(dy < 0)            +
                                               9*(abs(dx) < abs(dy)) +
                                              18*(abs(dx) > abs(dy)) ]


def winding(points, last_point = None):
    """Returns the winding number in octants around origin,
       or None if the polyline passes through origin."""
    total = 0

    octant_changes = (1, 2, 3, None, -3, -2, -1, 0, 1, 2, 3, None, -3, -2, -1)

    # If the last point is provided as a separate parameter, then points can be a generator.
    if last_point is None:
        next_point = points[-1]
    else:
        next_point = last_point
    next_octant = octant(next_point[0], next_point[1])
    if next_octant is None:
        return None

    # Iterate over the points
    for p in points:
        curr_point = next_point
        curr_octant = next_octant

        next_point = p
        next_octant = octant(next_point[0], next_point[1])
        if next_octant is None:
            return None

        # If no change in octants, nothing to do.
        if curr_octant == next_octant:
            continue

        # Limit/wrap octant change to [-4, +3] inclusive.
        change = octant_changes[7 + next_octant - curr_octant]
        if change is None:
            # Line segment makes almost a half a turn around origin.  Which way?
            delta = curr_point[0]*next_point[1] - curr_point[1]*next_point[0]
            if delta > 0:
                total += 4
            elif delta < 0:
                total -= 4
            else:
                return None
        else:
            total += change

    return total


def rotation_matrix(to_z_axis):
    # Error "loop": breaking out of this returns the identity matrix.
    while True:
        # Norm/length
        d = sqrt(float(to_z_axis[0])**2 + float(to_z_axis[1])**2 + float(to_z_axis[2])**2)
        if d <= 0.0:
            break

        # Unit vector to be rotated to (0,0,1):
        old_z = ( float(to_z_axis[0]) / d,
                  float(to_z_axis[1]) / d,
                  float(to_z_axis[2]) / d )

        # Rotation axis is (old_z[1], -old_z[0], 0.0).
        d = sqrt(old_z[0]*old_z[0] + old_z[1]*old_z[1])
        if d <= 0.0:
            break

        # Rotation axis unit vector x and y components
        x =  old_z[1] / d
        y = -old_z[0] / d

        # Cosine, sine, and 1 - cosine of the rotation angle:
        c = old_z[2]
        s = sqrt(1 - c*c)
        u = 1.0 - c

        # Rotation matrix based on unit axis and angle, with z=0:
        return (( c + x*x*u,     x*y*u,  y*s ),
                (     y*x*u, c + y*y*u, -x*s ),
                (      -y*s,       x*s,    c ))

    # Otherwise, return the unit matrix.
    return ((1,0,0),(0,1,0),(0,0,1))


def transform_3d_2d(point, relative, R):
    x = point[0] - relative[0]
    y = point[1] - relative[1]
    z = point[2] - relative[2]
    return ( R[0][0]*x + R[0][1]*y + R[0][2]*z,
             R[1][0]*x + R[1][1]*y + R[1][2]*z )

def winding3d(points, axis_point=(0.0, 0.0, 0.0), axis_direction=(0.0, 0.0, 1.0)):
    """Returns the winding number in octants around the specified point and axis,
       or None if the polyline passes through the axis."""
    octant_changes = (1, 2, 3, None, -3, -2, -1, 0, 1, 2, 3, None, -3, -2, -1)
    total = 0

    R = rotation_matrix(axis_direction)

    next_point = transform_3d_2d(points[-1], axis_point, R)
    next_octant = octant(next_point[0], next_point[1])
    if next_octant is None:
        return None

    for p in points:
        curr_point = next_point
        curr_octant = next_octant

        next_point = transform_3d_2d(p, axis_point, R)
        next_octant = octant(next_point[0], next_point[1])
        if next_octant is None:
            return None

        if curr_octant == next_octant:
            continue

        change = octant_changes[7 + next_octant - curr_octant]
        if change is None:
            delta = curr_point[0]*next_point[1] - curr_point[1]*next_point[0]
            if delta > 0:
                total += 4
            elif delta < 0:
                total -= 4
            else:
                return None
        else:
            total += change

    return total

if __name__ == '__main__':
    from sys import argv, exit, stderr

    if len(argv) < 6:
        stderr.write("\n")
        stderr.write("Usage: %s [ -h | --help ]\n" % argv[0])
        stderr.write("       %s Xa,Ya,Za Xp,Yp,Zp X0,Y0,Z0 X1,Y1,Z1 ... XN,YN,ZN\n" % argv[0])
        stderr.write("\n")
        stderr.write("Where:\n");
        stderr.write("       Xa,Ya,Za   is the axis direction around which we count winding,\n")
        stderr.write("       Xp,Yp,Zp   is a point on the axis, and\n")
        stderr.write("       Xi,Yi,Zi   are the points on the polyline.\n")
        stderr.write("\n")
        exit(0)

    points = []
    for p in argv[1:]:
        triplet = p.replace(',', ' ').replace(':', ' ').replace('/', ' ').strip().split()
        if len(triplet) != 3:
            stderr.write("%s: Not a 3D point.\n" % p)
            exit(0)
        else:
            points.append(tuple(map(float, triplet)))

    points = list(points)
    axis_direction = points.pop(0)
    axis_point = points.pop(0)

    w = winding3d(points, axis_point, axis_direction)
    if w is None:
        print("Polyline passes through the axis.")
    elif w == 0:
        print("The axis is outside the polyline.")
    elif w == 8 or w == -8:
        print("The axis is inside the polyline.")
    else:
        print("The winding number is %d." % w)

You can run it for example via

python3 -B example.py  0,0,1 0,0,0  2,1,0 1,2,0 -1,2,0 -2,1,0 -2,-1,0 -1,-2,0 1,-2,0 2,-1,0

which will output

The axis is inside the polyline.

because the Z axis at origin is inside the octagon. You can also import it, in which case it exports four functions:

• winding(points [, last_point):

  Returns the winding number (in octants) for the 2D polygon with vertices in points. Normally, points is a list or a tuple (of tuples), but if the last_point is specified (as the last point in the sequence), then point can be a generator/iterator.

• winding3d(points, axis_point, axis_direction):

  Returns the winding number (in octants) for the 3D polygon with vertices in points, with respect to the axis specified by direction (axis_direction) and a point on the axis (axis_point). The direction does not need to be an unit vector.

• octant(dx, dy)

  Returns the octant as specified in this answer, from $0$ to $7$, inclusive, or None if dx==0 and dy==0 (or if there is a bug in the implementation, although light testing indicates it should be correct).

• rotation_matrix(to_z_axis)

  Returns a tuple of tuples describing a 3×3 rotation matrix that rotates to_z_axis parallel to the (positive) $z$ axis.

Answer 2

Interesting problem! It can be done with some easy steps:

1. Pick a random point A on any edge (except vertex).
2. Find the equation of half-plane P formed by the red line and point A.
3. Find all intersections (except vertex) of edges and half-plane P.
4. If the number of intersections is odd, then the polyline is around the line.

P/s: check the better solution in @Example's comment