Saturday, 24 August 2013

Segment Intersection Demo

An addendum to my previous post in the quest to find a solution to polygon diameter, antipodal distance and/or polygon fetch for concave polygons.  This ditty checks to see if line segments intersect...may come in handy...read the header etc.

'''
 SegmentIntersectionDemo.py

 Author
  Dan Patterson
 
obidangis@gmail.com
  Dan_Patterson@carleton.ca
Purpose:
  Determines whether two line segments intersect and returns it if it does
  Detects parallel lines.  This differs from line intersections, whereby the
  point of intersection need not be on either line
 
Background:
 
http://gis.stackexchange.com/questions/69025/antipodal-distance-or-polygon \
  -fetch-or-polygon-diameter-for-concave-polygons
'''
def SegmentIntersection(A,B,C,D):
  #[read the point objects, A,B,C,D which are lists of X,Y coordinates
  x0,y0 = A; x1, y1 = B; x2,y2 = C; x3,y3 = D
  # calculate differences
  Ax = x1-x0;  Ay = y1-y0
  Bx = x3-x2;  By = y3-y2
  Cx = x0-x2;  Cy = y0-y2

  # calculate the lengths of the two lines
  len1 = math.sqrt(Ax*Ax + Ay*Ay)
  len2 = math.sqrt(Bx*Bx + By*By)

  # calculate angle between the two lines
  dot = (Ax*Bx + Ay*By) #dot product
  deg = dot/(len1*len2)

  #if abs(angle) == 1:  #lines are parallel, no intersection is possible
  if(abs(deg) == 1):
    return None

  #find the intersection point (pt) between two lines
  pt = [0,0]
  div = By*Ax - Bx*Ay
  ua = (Bx*Cy - By*Cx)/div
  ub = (Ax*Cy - Ay*Cx)/div
  pt[0] = x0 + ua*Ax
  pt[1] = y0 + ua*Ay

  # calculate the combined length of the two segments
  # between pt- pnt0 and pt-pnt1
  Ax = pt[0] - x0
  Bx = pt[0] - x1
  Ay = pt[1] - y0
  By = pt[1] - y1
  segmentLen1 = math.sqrt(Ax*Ax + Ay*Ay) + math.sqrt(Bx*Bx + By*By)
  #print "segmentLen1 ", segmentLen1
  # calculate the combined length of the two segments
  # between Pt-p2 and Pt-p3
  Ax = pt[0] - x2
  Bx = pt[0] - x3
  Ay = pt[1] - y2
  By = pt[1] - y3
  segmentLen2 = math.sqrt(Ax*Ax + Ay*Ay) + math.sqrt(Bx*Bx + By*By)

  '''
  if the lengths of both sets of segments are the same as
     the lengths of the two lines the point is actually
   on the line segment.
   if the point isn't on the line, return None, else the intersection
  '''
  if(abs(len1-segmentLen1) > 0.0 or abs(len2-segmentLen2) > 0.0):
    return None
  else:
    return pt

if __name__ == "__main__":
  import math
  print "\nSegment Intersection Demo...Tests to see if two line segments intersect"
  dataList = [[[0.0,0.0],[1.0,1.0], [0.0,1.0],[1.0,0.0]],
              [[0.0,0.0], [0.0,1.0],[0.0,1.0],[1.0,1.0]],
              [[0.0,0.0], [0.0,1.0],[1.0,1.0],[0.0,1.0]],
              [[0.0,0.0],[0.25,0.75],[0.0,1.0],[1.0,1.0]]
              ]
  print dataList[0]
  for pnts in dataList:
    A,B,C,D = pnts
    pt = SegmentIntersection(A,B,C,D)
    print "\nintersection point = ", pt, " for pnts: ",A,B,C,D

No comments:

Post a Comment