Algorithms for calculating vertex normals in 2D polygon

Question

I have the points of a 2D polygon in the order in which they are connected, and I need to calculate the vertex normals and move them along it. An example is shown.

I would like to know how to calculate them efficiently and which structures I could use. I don't know if I have to find the orientation of the polygon, and what to do in general. As I said the only thing I have to start with is the list of points in the order in which they are connected. Say 5 vertices, I have a list that goes: x1, y1, x2, y2, x3, y3, x4, y4, x5, y5.

Where should I start?

It doesn't matter if the polygon becomes complex or self crossing, nothing matters just that to calculate the normals at the vertices and move them along these normals in the correct direction, all outwards or all inwards the polygon.

Answer 1

If I understand the question correctly we have the following:

• Given a list of numbers (I'm going to presume we're working with doubles or decimals) we need to transform this into a sequence of 2D vectors which describe a path along the perimeter of a polygon either clockwise or anticlockwise (which we can know ahead of time)

• Next we calculate the vector between each pair of points, including the first and last. We're only interested in the direction of this vector rather than the size, so we may as well make these unit vectors by dividing them by their magnitude.

• So now we have two same-length sequences of vectors, one describing nodes, the other the direction along each edge. If we take these tangential vectors we can rotate them by a quarter-turn in order to get vectors in the normal direction to each edge (whether you go clockwise or anticlockwise depends on which way around the shape your points go), and the sum + normalise each pair, again including the first and last.

• What we have now is a sequence of unit vectors at a normal to each vertex, which we can add or subtract in order to grow or shrink our shape.

Looking at the operations we need to perform leads us towards a data-structure. I would propose that we will obviously want something to represent a 2d vector, and it makes sense for this to be immutable for reasons I'll describe later. I would be tempted to create another data structure for our shape, constructed with a doubly-linked-list of vectors, (doubly-linked for reasons which should become clear later) which exposes a same-length sequence of pairs of vectors to represent a point and its normal. These are standard enough data-structures that I'm not going to provide pseudo-code to implement them.

This constructor can then hold a dictionary/hashmap of vector pair to vector - this is constructed from our rotated unit tangential vectors keyed by the two vertices which form the described edge. Having vectors immutable allows for good behaviour using their references/pointers as keys allowing for efficient lookup later.

We can then for each point in our input create a pair of this point, and the sum of the dictionary entries keyed by this and its two neighbours (which is why doubly-linked), which is a vector along the normal, which we then divide by its magnitude to make this a unit vector too.

I have made some assumptions here, including that you're using an OO language, but similar patterns can be employed just as effectively in other paradigms. The other assumptions I made:

• The points given follow either clockwise or anticlockwise - if they are not in some order this is rather insoluble as you can't know what shape is being described,

• That you know whether it's clockwise or anticlockwise up front - if not this can be deduced by summing the angles to the left and right as you go round (not forgetting to join the first and last). If the angles to the right are greater, you are traveling anticlockwise, if it's left, clockwise.

Code example of the above waffle:

""" 
For the purposes of pseudo-code I'm using python for its clarity and I have
gone for a functional style for its brevity. 

I'm assuming you have a `Vector2D` class which supports vector addition + 
scalar division and has a property `normalised` which returns the result of 
dividing the vector by its magnitude.

I'm also assuming you're using an implementation of a doubly linked list that 
is something like https://dbader.org/blog/python-linked-list

`normal_of` will be some matrix multiplication of a vector, to perform the 
quarter rotation, but which direction this is depends on the input.
"""

@dataclass
class VertexInfo:
    vertex: Vector2D
    normal: Vector2D

def unit_direction_between(from_vertex: Vector2D, to_vertex: Vector2D) -> Vector2D:
    return (to_vertex - from_vertex).normalised

def unit_average_vectors(vec1: Vector2D, vec2: Vector2D) -> Vector2D:
    return (vec1 + vec2).normalised

def iterate_pairs_from(vertices: DoublyLinkedList[Vector2D]) -> Iterator[Tuple[Vector2D, Vector2D]]:
    """ This covers the fact that we need to associate the first and last """
    node = first = vertices.head
    while node:
        following = node.next
        yield node.value, (following or first).value
        node = following

def shape_from_vertices(vertices: DoublyLinkedList[Vector2D]) -> Iterator[VertexInfo]:
    edges = {
        (node, following): normal_of(unit_direction_between(node, next))
        for (node, following) in iterate_pairs_from(vertices)
    }

    # A bit of a dance to get the last entry in the list, this would be unnecessary with 
    # a standard python list (which is actually a dynamic array) or tuple as you can
    # simply index to collection[-1] - this could be implemented on your doubly 
    # linked list if you track the last node on your list object.
    node = first = vertices.head
    while node:
        node, last = node.next, node

    # This is the actual business logic
    node = first
    while node:
        previous = (node.previous or last).value
        current = node.value
        following = (node.next or first).value

        edge_normal_behind = edges[previous, current]
        edge_normal_ahead = edges[current, following]
        vertex_normal = unit_average_vectors(edge_normal_behind, edge_normal_ahead)

        yield VertexInfo(current, vertex_normal)

        node = node.next

There are some parts of this which could be simplified - you don't really need a dictionary or a doubly-linked-list but they make the code a lot more obvious without any significant performance penalties, and relying on well-known data structures makes it easier to test your implementation.

Answer 2

First, obtain the normals for each edge of the polygon. The normal for an edge is given by the normalized cross product of the edge vector (p2 - p1) with the 2D plane normal (a unit vector pointing in the direction of the Z-axis). The cross product of a 2D vector with the positive Z-axis is given by (-y, x). This will give you outward pointing normals for vertices in the counter-clockwise winding order, and inward pointing normals for a clockwise vertex winding order. Let the normals be stored in a list of n (number of vertices) elements, where the normal of the edge formed by vertices i and (i + 1) % n is stored at the index i.

For each vertex i, its neighboring edge normals have the indices i and (i - 1) % n. The simplest way to combine these edge normals into a vertex normal would be a normalized sum ((normals[i] + normals[(i - 1) % n]).normalized).

In pseudocode:

normals = [(0, 0)] * n

for i in 0..n {
    j = (i + 1) % n
    edge_vec = vertices[j] - vertices[i]
    normals[i] = (-edge_vec.y, edge_vec.x).normalized
}

vertex_normals = [(0, 0)] * n


for i in 0..n {
    j = (n + i - 1) % n
    vertex_normals[i] = (normals[i] + normals[j]).normalized
}

For the vertex normal calculation algorithm, more "intuitive" vertex normals may be calculated using a weighted edge normal sum instead, where each edge normal is weighted by the edge length (normals[i] * edges[i].length + normals[j] * edges[j].length). Since edge length is equal to the magnitude of the edge vector which is the same as the magnitude of the un-normalized edge normal from our cross product, we can skip normalization of the normals altogether, resulting in a (theoretically) faster algorithm.

weighted_normals = [(0, 0)] * n

for i in 0..n {
    j = (i + 1) % n
    edge_vec = vertices[j] - vertices[i]
    weighted_normals[i] = (-edge_vec.y, edge_vec.x) // == normals[i] * edge_vec.length
}

vertex_normals = [(0, 0)] * n


for i in 0..n {
    j = (n + i - 1) % n
    vertex_normals[i] = (weighted_normals[i] + weighted_normals[j]).normalized
}

Your animated example seems to be using a weighted sum too, as the vertex normal of the bottom-most concave angled vertex stays the same even though the edge normal to its left changes in direction while the edge normal to the right stays almost the same.

Answer 3

What you call the vertex normals are not normals. They are the bissectors of the normals of the two adjoining sides.

To find their direction, normalize the normals to unit length and take the vector sum. You can normalize once again. Make sure to take the normals on the same side of the edges.

Answer 4

First compute the edge normals as follows. Treat each edge e as a vector from one endpoint to the other (subtracting one endpoint from another). Then rotate edge e 90 degrees by multiplying by the rotation matrix. $$ R = \left( \begin{array}{rr} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array} \right) = \left( \begin{array}{rr} 0 & 1 \\ -1 & 0 \end{array} \right) $$ Now normalize this edge normal to have unit length.

Then, as @DMGregory suggested, average two edge normals a and b by forming a + b and normalizing to unit length.