Python & Algorithm: How to do simple geometry shape match?

Question

Given a set of points (with order), I want to know if its shape is within certain types. The types are:

rectangle = [(0,0),(0,1),(1,1),(1,0)]
hexagon = [(0,0),(0,1),(1,2),(2,1),(2,0),(1,-1)]
l_shape = [(0,0),(0,3),(1,3),(1,1),(3,1),(3,0)]
concave = [(0,0),(0,3),(1,3),(1,1),(2,1),(2,3),(3,3),(3,0)]
cross = [(0,0),(0,-1),(1,-1),(1,0),(2,0),(2,1),(1,1),(1,2),(0,2),(0,1),(-1,1),(-1,0)]

For example, give roratated_rectangle = [(0,0),(1,1),(0,2),(-1,1)] We'll know it belongs to the rectangle above.

is similar to

Notice:

1. rotation and edge with different length are considered similar.
2. Input points are ordered. (so they can be draw by path in polygon module)

How can I do it? Is there any algorithm for this?

What I'm thinking:

Maybe we can reconstruct lines from given points. And from lines, we can get angles of a shape. By comparing the angle series (both clockwise and counterclockwise), we can determine if input points are similar to the types given above.

Answer 1

Your thinking is basically the right idea. You want to compare the sequence of angles in your test shape to the sequence of angles in a predefined shape (for each of your predefined shapes). Since the first vertex of your test shape may not correspond to the first vertex of the matching predefined shape, we need to allow that the test shape's sequence of angles may be rotated relative to the predefined shape's sequence. (That is, your test shape's sequence might be a,b,c,d but your predefined shape is c,d,a,b.) Also, the test shape's sequence may be reversed, in which case the angles are also negated relative to the predefined shape's angles. (That is, a,b,c,d vs. -d,-c,-b,-a or equivalently 2π-d,2π-c,2π-b,2π-a.)

We could try to choose a canonical rotation for the sequence of angles. For example, we could find the lexicographically minimal rotation. (For example, l_shape's sequence as given is 3π/2,3π/2,π/2,3π/2,3π/2,3π/2. The lexicographically minimal rotation puts π/2 first: π/2,3π/2,3π/2,3π/2,3π/2,3π/2.)

However, I think there's a risk that floating-point rounding might result in us choosing different canonical rotations for the test shape vs. the predefined shape. So instead we'll just check all the rotations.

First, a function that returns the sequence of angles for a shape:

import math

def anglesForPoints(points):
    def vector(tail, head):
        return tuple(h - t for h, t in zip(head, tail))

    points = points[:] + points[0:2]
    angles = []
    for p0, p1, p2 in zip(points, points[1:], points[2:]):
        v0 = vector(tail=p0, head=p1)
        a0 = math.atan2(v0[1], v0[0])
        v1 = vector(tail=p1, head=p2)
        a1 = math.atan2(v1[1], v1[0])
        angle = a1 - a0
        if angle < 0:
            angle += 2 * math.pi
        angles.append(angle)
    return angles

(Note that computing the cosines using the dot product would not be sufficient, because we need a signed angle, but cos(a) == cos(-a).)

Next, a generator that produces all rotations of a list:

def allRotationsOfList(items):
    for i in xrange(len(items)):
        yield items[i:] + items[:i]

Finally, to determine if two shapes match:

def shapesMatch(shape0, shape1):
    if len(shape0) != len(shape1):
        return False

    def closeEnough(a0, a1):
        return abs(a0 - a1) < 0.000001

    angles0 = anglesForPoints(shape0)
    reversedAngles0 = list(2 * math.pi - a for a in reversed(angles0))
    angles1 = anglesForPoints(shape1)
    for rotatedAngles1 in allRotationsOfList(angles1):
        if all(closeEnough(a0, a1) for a0, a1 in zip(angles0, rotatedAngles1)):
            return True
        if all(closeEnough(a0, a1) for a0, a1 in zip(reversedAngles0, rotatedAngles1)):
            return True
    return False

(Note that we need to use a fuzzy comparison because of floating point rounding errors. Since we know the angles will always be in the small fixed range 0 … 2π, we can use an absolute error limit.)

>>> shapesMatch([(0,0),(1,1),(0,2),(-1,1)], rectangle)
True
>>> shapesMatch([(0,0),(1,1),(0,2),(-1,1)], l_shape)
False
>>> shapesMatch([(0,0), (1,0), (1,1), (2,1), (2,2), (0,2)], l_shape)
True

If you want to compare the test shape against all the predefined shapes, you may want to compute the test shape's angle sequence just once. If you're going to test many shapes against the predefined shapes, you may want to precompute the predefined shapes' sequences just once. I leave those optimizations as an exercise for the reader.

Answer 2

Your shapes may not only be rotated, they are also translated and may even be scaled. The orientation of the nodes may also differ. Fo example your original square has an edge length of 1.0 and is defined anticlockwise, whereas your diamond shape has an edge length of 1.414 and is defined clockwise.

You need to find a good reference for comparison. The following should work:

• Find the centre of gravity C for each shape.
• Determine radial coordinates (r, φ) for all nodes, where the origin of the radial coordinate system is the centre of gravity C.
• Normalise the radii for each shape so that the maximum value of r in a shape is 1.0
• Ensure that the nodes are oriented anticlockwise, i.e. with increasing angles φ.

Now you have two lists of n radial coordinates. (The cases where the number of nodes in a shape don't match or where there are fewer than three nodes should have been ruled out already.)

Evaluate all n configurations of offsets, where you leave the first array as it is and shift the second. For a four-element array, you compare:

{a1, a2, a3, a4} <=> {b1, b2, b3, b4}
{a1, a2, a3, a4} <=> {b2, b3, b4, b1}
{a1, a2, a3, a4} <=> {b3, b4, b1, b2}
{a1, a2, a3, a4} <=> {b4, b1, b2, b3}

The radial coordinates are floating-point numbers . When you compare values, you should allow for some latitude to cater for inaccuracies introduced by the floating-point math. Because the numbers 0 ≤ r ≤ 1 and −π ≤ φ ≤ π are roughly in the same range, you can used a fixed epsilon for that.

Radii are compared by their normalised value. Angles are compared by their difference to the angle of the previous point in the list. When that difference is negative, we have wrapped around the 360° border and have to adjust it. (We must enforce positive angles differences, because the shape we compare to might not be equally rotated and thus may not have the wrap-around gap.) The angle is allowed to go both forwards and backwards, but must eventually come full circle.

The code has to check n configurations and test all n nodes for each. In practice, mismatches will be found early, so that the code should have okay performance. If you are going to compare a lot of shapes, it might be worth creating the normalised, anticlockwise radial representation for all shapes beforehand.

Anyway, here goes:

def radial(x, y, cx = 0.0, cy = 0.0):
    """Return radial coordinates from Cartesian ones"""

    x -= cx
    y -= cy

    return (math.sqrt(x*x + y*y), math.atan2(y, x))



def anticlockwise(a):
    """Reverse direction when a is clockwise"""

    phi0 = a[-1]
    pos = 0
    neg = 0

    for r, phi in a:
        if phi > phi0:
            pos += 1
        else:
            neg += 1

        phi0 = phi

    if neg > pos:
        a.reverse()



def similar_r(ar, br, eps = 0.001):
    """test two sets of radial coords for similarity"""

    _, aprev = ar[-1]
    _, bprev = br[-1]

    for aa, bb in zip(ar, br):
        # compare radii
        if abs(aa[0] - bb[0]) > eps:
            return False

        # compare angles
        da = aa[1] - aprev
        db = bb[1] - bprev

        if da < 0: da += 2 * math.pi
        if db < 0: db += 2 * math.pi

        if abs(da - db) > eps:
            return False

        aprev = aa[1]
        bprev = bb[1]

    return True



def similar(a, b):
    """Determine whether two shapes are similar"""

    # Only consider shapes with same number of points
    if len(a) != len(b) or len(a) < 3:
        return False        

    # find centre of gravity
    ax, ay = [1.0 * sum(x) / len(x) for x in zip(*a)]
    bx, by = [1.0 * sum(x) / len(x) for x in zip(*b)]

    # convert Cartesian coords into radial coords
    ar = [radial(x, y, ax, ay) for x, y in a]
    br = [radial(x, y, bx, by) for x, y in b]

    # find maximum radius
    amax = max([r for r, phi in ar])
    bmax = max([r for r, phi in br])

    # and normalise the coordinates with it
    ar = [(r / amax, phi) for r, phi in ar]
    br = [(r / bmax, phi) for r, phi in br]

    # ensure both shapes are anticlockwise
    anticlockwise(ar)
    anticlockwise(br)

    # now match radius and angle difference in n cionfigurations
    n = len(a)
    while n:
        if similar_r(ar, br):
            return True                

        br.append(br.pop(0))      # rotate br by one
        n -= 1

    return False

Edit: While this solution works, it is overcomplicated. Rob's answer is similar in essence, but uses a straightforward metric: the internal angles between the edges, which takes care of translation and scaling automatically.

Answer 3

Approaching this by linear algebra, each point will be a vector (x,y), which can be multiplied into a rotation matrix to get new coordinates (x1,y1) by the designated angle. There does not appear to be LaTeX support here so I can't write it out clearly, but in essence:

(cos(a) -sin(a);sin(a) cos(a))*(x y) = (x1 y1)

This results in coordinates x1, y1 that are rotated by angle "a".

EDIT: This is perhaps the underlying theory, but you can set up an algorithm to shift the shape point-by-point to (0,0), then compute the cosine of the angle between adjacent points to classify the object.