8

I am working on some class examples in Python implemented within ArcMap to calculate the antipodal distance within a polygon. This is fairly routine for convex polygons, however, for concave polygons, I wish to exclude solutions (formed by a ray connecting the boundary points), that aren't completely within the polygon and not on the polygon boundary or intersecting it. Have I interpreted the definition wrong or is this beast go by another name.

Consider these two polygons

pnts = [ [0,0], [0,1], [1,4],[3,5],[5,4],[4,1],[0,0] ] # a closed-loop convex

pnts = [ [0,0], [2,1], [1,4],[3,5],[5,4],[4,1],[0,0] ] # a closed-loop concave polygon

In my interpretation, point 0,0 should not have an antipodal distance associated with it since the vector connecting it with the other points is either self intersecting the polygon or is on the polygon boundary.

If anyone has any clarification on the definition or potential solutions, I would appreciate it.

A visual of the convex polygon and the desired lines (shown in red) is enclosed (example vectors from point 0 are only shown).

Interior Antipodal Distance example

In the convex example, the first point has no antipodal vectors, however, the second point does.

Concave Antipodal Example

EDIT I have had some success searching using "polygon fetch" or "polygon diameter" on the web, I suspect that this is what I am after.

  • 1
    Hi, Dan. What definition of "antipodal distance" are you using? One possibility would be the furthest point as measured by travel along the polygon's boundary, but that does not seem consistent with your description. Another definition is a furthest point where travel can occur anywhere inside or outside the polygon. Yet a third is the furthest point where travel is allowed only within the interior and boundary of the polygon. – whuber Aug 19 '13 at 14:28
  • 1
    @whuber, I was looking for a solution which only travelled within the polygon excluding line segments which form the polygon boundary. In the convex example I gave, movement from points p0 to p1, or p0 to p5 would not be permitted since they are part of the polygon edge, however, p0 to p2, p3, p4 would be. Hence, my concern that "antipodal" may not be the correct term. Note, I am only interested in single-part convex polygons with no holes at the moment. If I am stuck with edge segments in the solution, I can always remove them later. – user681 Aug 19 '13 at 15:00
  • 1
    There is a delicate issue here, Dan: although such segments might be ruled out, nevertheless they do tell you what the infimum of all possible distances will be (they merely prevent that infimum from actually being realized). Practical solutions would keep on the inside of those segments but remain infinitesimally close to them. Thus, for convex polygons the algorithm is simple: find a vertex furthest from the starting point (there can be many of them: imagine a semicircle and starting at the center of the original circle). – whuber Aug 19 '13 at 15:09
  • 1
    I still don't understand your definition, Dan, because there is no "longest path" within any polygon: you can snake around to make arbitrarily long paths. Possibly what you intend is the following: define the distance between points P and Q in a (connected) polygon to be the infimum of the lengths of all paths from P to Q lying wholly within the polygon. Then a plausible "antipode" for a compact connected polygon P would be any point Q at the maximum distance from P. (When P is a vertex of a convex polygon its antipodes once again are vertices at the maximum Euclidean distance from P.) – whuber Aug 20 '13 at 19:39
  • 2
    The furthest point is rigorously characterized using the "plausible" definition in my preceding comment. Note that in finding it you are allowed to assume that you may travel along edges. In your second figure, E is antipodal to A and B; A is antipodal to C, D, and E; and D and A are both antipodal to F. Using the canoe analogy, where the interior of the polygon is a lake, a point P is antipodal to your starting point Q when in a canoe race from Q against an opponent who aims to reach P before you reach some point P', they have no advantage over you no matter where P' is. – whuber Aug 21 '13 at 19:08
4
+100

If I were writing an algorithm I would simply check if a line between two vertices on the polygon intersects any line that forms one of the edges. Here's my pseudo code:

  1. Identify all the vertices, store in a list
  2. Identify all the edges, store in a list (as pairs of vertices from 1, perhaps)
  3. for each vertex, get the distance to all other vertices, except:
    1. exclude neighbouring vertices (those that share a pairing with this vertex in 2, perhaps)
    2. exclude any line that intersects any line in 2. using something from here.
  4. store all the valide distances with reference to the vertices in 1.

  5. do what ever you want with the results, write new lines out, store the longest one for each polygon ...

Now, I'm not sure if this is what you're after, but you certainly can do the above in ArcPy.

EDIT: Code for step 2.2:

E = B-A = ( Bx-Ax, By-Ay )
F = D-C = ( Dx-Cx, Dy-Cy ) 
P = ( -Ey, Ex )
h = ( (A-C) * P ) / ( F * P )

If h is between 0 and 1, the lines intersect, otherwise they don't. If F*P is zero, of course you cannot make the calculation, but in this case the lines are parallel and therefore only intersect in the obvious cases. If h is 1, then the lines end at the same point. Handle this as you will! (I'd say they intersect, it makes i easier.)

Another example for step 2.2 from here: http://en.wikipedia.org/wiki/Line-line_intersection enter image description here

First check the denominator doesn't equal 0, which means that the lines are parallel.

Then check that the coordinates found above are not outside the bounding box of the either line.

More reading: http://compgeom.cs.uiuc.edu/~jeffe/teaching/373/notes/x06-sweepline.pdf

  • The code on my blog does everything except 3.2 and refers the result back to the calling program for further calculations which work well for convex polygons. I would like some references if possible and whether the winding number would be efficient to determine line intersections or I should go some other route. – user681 Aug 22 '13 at 0:42
  • 2
    I added an example of the intersection calculation from the article I linked to. I think that I wouldn't worry about efficiency until it was an issue! Get something working, then fix it if it's not good enough! – Alex Leith Aug 22 '13 at 2:05
  • Thanks Alex, I will check it out...efficiency isn't an issue since this is only an example not to be run on thousands of polygons – user681 Aug 22 '13 at 2:15
  • Although it's hard to tell what this answer describes, it seems to be a check whether a polygon is convex. It's particularly confusing that it appears to use "line" in place of "line segment." The code given does not seem to be a valid check for intersection of line segments. – whuber Aug 22 '13 at 15:49
  • Hi whuber, this does check line segments. The second example I added perhaps makes it clearer, in that the math works out the intersectino point for infinite lines and you then need to check whether that intersection point is inside the bounding box of one of the lines. Like all vector math, there's surely a library that does this, but it's not that complex, and I think necessary for what OP wants to do. – Alex Leith Aug 23 '13 at 0:05
3

I would be tempted to do this using angels, almost like line of sight. If while iterating the vertices in the shape the angles between the origin vertex and destination vertex continue in a consistent direction, all the points are candidates for the antipodal. If an angle switches direction, then that point is hidden by or hides the previous point. If it is hidden by the previous point, the point needs to be skipped. If it hides the previous point, the previous point(s) need to be removed from the candidate list.

  1. Create a PolygonCandidates list
  2. For each vertex (point k)
    1. Create new list for candidates (Point, Angle)
    2. Add the current vertex to candidates list (point k)
    3. Iterate clockwise around polygon, for each remaining vertex (point i)
      1. If the angle to the current point (from point k to point i) continues in a clockwise direction 1.add the point
      2. If the angle to the current point continues in a counter-clockwise direction
      3. If the previous two candidate points, plus the current point form a right hand turn.
      4. Remove the last point in the list until the current angle and last candidate list angle is in a counterclockwise direction.
      5. Add the current point to the candidates list
    4. Add all but the first two, and last candidate points to a PolygonCandidates list
  3. Find the furthest point in the PolygonCandidates list.

I'm not sure what to do with cases where the origin, and two other vertices all fall along the same line. In that case, the angle would be the same. If you had a polygon with holes, you could find the min/max angle of each hole, and remove any candidate point that lies within that range.

The main advantage to this approach would be that you don't have to test for line intersection between the current line segment and all the polygon edges.

This works...I think. I have updated the pseudo code above and the python in order to make it easier to read.


This should be the last edit. The example below should find the largest anitpole for a given geometry. I altered the scrip so that is uses Points and Vectors, to try and make it easier to read.

import math
from collections import namedtuple


Point = namedtuple("Point", "position x y")
Vector = namedtuple("Vector", "source dest angle")

def isClockwise(angle1, angle2):
    diff = angle2 - angle1
    #print("         angle1:%s angle2:%s diff: %s" % (angle1, angle2, diff))
    if(diff > math.pi/2):
        diff = diff - math.pi/2
    elif (diff < -math.pi/2):
        diff = diff + math.pi/2
    #print("         diff:%s" % (diff)) 
    if(diff > 0):
        return False
    return True

def getAngle(origin, point):
    return math.atan2(point.y - origin.y, point.x-origin.x)

#returns a list of candidate vertcies.  This will include the first, second, and second to last points 
#the first and last points in the polygon must be the same
#k is the starting position, only vertices after this position will be evaluated
def getCandidates (k, polygon):

    origin = polygon[k]
    candidates = [Vector(k,k,0)]
    prevAngle = 0;
    currentAngle = 0;
    for i in range(k + 1, len(polygon) - 1):

        current = polygon[i]
        #print("vertex i:%s x:%s y:%s  " % (i, current.x, current.y))

        if(i == k+1):
            prevAngle = getAngle(origin, current)
            candidates.append(Vector(k,i,prevAngle))
        else:   
            currentAngle = getAngle(origin, current)
            #print("     prevAngle:%s currentAngle:%s  " % (prevAngle, currentAngle))
            if isClockwise(prevAngle, currentAngle):
                #print("     append")
                candidates.append(Vector(k,i,currentAngle))
                prevAngle = currentAngle
            else:
                #look at the angle between current, candidate-1 and candidate-2
                if(i >= 2):
                    lastCandinate = polygon[candidates[len(candidates) - 1].dest]
                    secondLastCandidate = polygon[candidates[len(candidates) - 2].dest]
                    isleft = ((lastCandinate.x - secondLastCandidate.x)*(current.y - secondLastCandidate.y) - (lastCandinate.y - secondLastCandidate.y)*(current.x - secondLastCandidate.x)) > 0
                    #print("     test for what side of polygon %s" % (isleft))
                    if(i-k >= 2 and not isleft):
                        while isClockwise(currentAngle, candidates[len(candidates) - 1].angle):
                            #print("     remove %s" % (len(candidates) - 1))
                            candidates.pop()
                        #print("     append (after remove)")
                        candidates.append(Vector(k,i,currentAngle))
                        prevAngle = currentAngle

        #for i in range(len(candidates)):
        #   print("candidate i:%s x:%s y:%s a:%s " % (candidates[i][0], candidates[i][1], candidates[i][2], candidates[i][3]))

    return candidates

def calcDistance(point1, point2):
    return math.sqrt(math.pow(point2.x - point1.x, 2) + math.pow(point2.y - point1.y, 2))

def findMaxDistance(polygon, candidates):
    #ignore the first 2 and last result
    maxDistance = 0
    maxVector = Vector(0,0,0);
    for i in range(len(candidates)):
        currentDistance = calcDistance(polygon[candidates[i].source], polygon[candidates[i].dest])
        if(currentDistance > maxDistance):
            maxDistance = currentDistance
            maxVector = candidates[i];
    if(maxDistance > 0):
        print ("The Antipodal distance is %s from %s to %s" % (maxDistance, polygon[candidates[i].source], polygon[candidates[i].dest]))
    else:
        print ("There is no Antipodal distance")

def getAntipodalDist(polygon):
    polygonCandidates = []
    for j in range(0, len(polygon) - 1):
        candidates = getCandidates(j, polygon)
        for i in range(2, len(candidates) - 1):
            #print("candidate i:%s->%s x:%s y:%s  " % (candidates[i].source, candidates[i].dest, candidates[i].x, candidates[i].y))
            polygonCandidates.append(candidates[i])

    for i in range(len(polygonCandidates)):
        print("candidate i:%s->%s" % (polygonCandidates[i].source, polygonCandidates[i].dest))
    findMaxDistance(polygon, polygonCandidates)


getAntipodalDist([Point(0,0,0),Point(1,-2,0),Point(2,-2,3),Point(3,2,2),Point(4,-1,1),Point(5,4,0),Point(6,0,0)])
getAntipodalDist([Point(0,0,0),Point(1,2,1),Point(2,1,4),Point(3,3,5),Point(4,5,4),Point(5,4,1),Point(6,0,0)])
getAntipodalDist([Point(0,0,0),Point(1,1,1),Point(2,2,1),Point(3,1,4),Point(4,3,5),Point(5,5,4),Point(6,4,1),Point(7,0,0)])
getAntipodalDist([Point(0,0,0),Point(1,-1,3),Point(2,1,4),Point(3,3,3),Point(4,2,0),Point(5,-2,-1),Point(6,0,0)])
  • Does this algorithm really work? Could you perhaps illustrate it with a simple example, such as the hexagon ((0,0),(-2,0),(-2,3),(2,2),(-1,1),(4,0),(0,0))? What happens when you start off at (0,0)? And what does this algorithm purport to do? For instance, it does not find the longest line segment in the polygon (of length 1.2*sqrt(26)). – whuber Aug 22 '13 at 20:47
  • Thanks for the comment Travis, however, this won't work in all cases (see concave hull example), from isRightTurn(A,B,C) would be false, and AC would not be a candidate segment. If B was further north, it could conceivably all one to for a segment AE, so I wouldn't want to rule out point A completely until all other points were checked. – user681 Aug 22 '13 at 20:48
  • @whuber, given that geometry, I don't see how the longest line segment is 1.2*sqrt(26). Unless I have totally missed what this question is about. Wouldn't it be sqrt(2), from either (0,0)->(-1,1) or (-2,0)->(-1,1). – travis Aug 23 '13 at 23:28
  • 1
    @DanPatterson, I may have missed what you are asking. My understanding was: whats the largest distance between a given vertex and any other vertex, that does not intersect the boundary of the polygon. I updated my script to find the polygon's maximum distance. – travis Aug 23 '13 at 23:34
  • Convex polygons don't seem to be an issue given the simplistic examples that one can find on the web and in texts. The polygon diameter for concave hulls seems to have some various interpretations and polygon fetch, I am now beginning to realize it another kettle of fish. In any event, non-intersecting whatever, is what I am after. My concern is my lack of clear definitions and examples with real-world examples. I can cover the convex ones, but the concave ones are proving problematic and beyond my mathematical/computational expertise as supported/highlighted by some of Bill's suggestions. – user681 Aug 23 '13 at 23:48
1

Perhaps consider triangulating the dataset. Which lines are common to polygons edges would be easy to establish and the remaining ones could be compared to find the longest? The question then is what triangulation algorithm you need.

It is only a hunch but I suspect (ironically) the "lowest quality" triangulation one can create must contain the line you are looking for e.g. Fig 1 in https://www.google.co.uk/url?sa=t&rct=j&q=&esrc=s&source=web&cd=6&ved=0CEoQFjAF&url=http%3A%2F%2Fhrcak.srce.hr%2Ffile%2F69457&ei=alIcUsb6HsLnswbfnYHoDw&usg=AFQjCNHIaykVRBAvv9hlaFJIBlfPLGHKtQ

  • My code on my blog does effectively than, it the terminology that I need clarified as well as what to do in the case of a concave hull. – user681 Aug 27 '13 at 23:51
  • A triangulation will inherently handle a concave hull (in as much as it will not create any triangle edges which cross a polygon boundary) – AnserGIS Aug 28 '13 at 18:30

Your Answer

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