Geometry Puzzler (part 2: intersecting segments)

[Mike Brenner]

Mike Fletcher > 2) Given two line segments, how can I find the
coordinates of their intersection, if there is one? (easy for right
angles, but I will often NOT have right angles.)

# This python software demonstrates how to compute the
# intersection of two line segments with no warranty.
# It is followed by a tiny test that shows how to call it. 
# It seems like all that kind of simple Algebra/Geometry
# should be available as part of the basic tools of the trade. 

# Most of it is just object oriented expression 
# that could be filtered out, 
# if this was the only thing that was being done 
# with the segments being intersected.

# It first computes the intersection of the lines on which the segment
lies, 
# if any, and then determines if that intersection lies on the first
segment, 
# thus there is no need to actually compute the second parameter.


import math

line_threshold=0.00001; # Closer than this, two lines are considered to
be the same line.

def tell(message):
  sys.stdout.softspace=0
  print(message);
         
def telln(message):
  print message;

def close(left, right, threshold):
  return abs(left-right)<=threshold;
  
class points2D:
  'Points === <x,y>'
  x_meters=0.0;
  y_meters=0.0;
  def __init__(point,x_meters,y_meters): # Point constructor,
    point.x_meters=float(x_meters)
    point.y_meters=float(y_meters)
  def x_of(point): # Extracts the x coordinate of a point
    return point.x_meters;
  def y_of(point):
    return point.y_meters; # Extracts the y coordinate of a point
  def __add__(left,right): # Adds two points
    assert(right.is_point2D());
    a=left.x_of()+right.x_of();
    b=left.y_of()+right.y_of();
    return points2D(a,b); 
  def __sub__(left,right): # Subtracts two points
    assert(right.is_point2D());
    a=left.x_of()-right.x_of();
    b=left.y_of()-right.y_of();
    return points2D(a,b); 
  def close(left,right,thresh): # Determines if left and right are close
to each other.
    assert(right.is_point2D());
    return close(left.x_of(),right.x_of(),thresh) and \
           close(left.y_of(),right.y_of(),thresh)
  def distance_squared (left, right): # Square of the distance between
two points.
    assert(right.is_point2D());
    return (right.x_of()-left.x_of())**2 + \
           (right.y_of()-left.y_of())**2
  def distance (left, right): # Distance between two points.
    assert(right.is_point2D())
    return math.sqrt(left.distance_squared(right));   
  def image(point): # Displayable image of a point.
    x,y=point.x_of(),point.y_of()
    return "["+`x`+", "+`y`+"]"
  def is_point2D(point): # Type test.
    return 1;
  def slope(left,right,threshold=0.00001): # Slope of a line.
    assert(right.is_point2D());
    denominator=right.x_of()-left.x_of();
    numerator  =right.y_of()-left.y_of();
    if close(denominator,0.0,threshold):
      telln("       slope left="+left.image())
      telln("            right="+right.image())
      raise slope_is_infinite_line_almost_vertical;
    return numerator/float(denominator);
  def times_scalar (point,scalar): # Scalar vector multiplication.
    a=scalar*point.x_of();
    b=scalar*point.y_of();
    return points2D(a,b);
  def y_intercept(point,slope): # The y_intercept of a nonvertical line.
    return point.y_of()-slope*point.x_of();
  
class lines:
  'Lines === <left_point,right_point>'
  left_point=points2D(0.0,0.0);
  right_point=points2D(0.0,0.0);
  def __init__(line,left_point,right_point): # Line constructor.
    assert(left_point.is_point2D())
    assert(right_point.is_point2D())
    if left_point.distance_squared(right_point)< line_threshold:
      raise lines_require_two_distinct_points
    line.left_point=left_point
    line.right_point=right_point
  def left_of(line): # Extracts the left point of a line.
    return line.left_point
  def right_of(line): # Extracts the right point of a line.
    return line.right_point
  def image(line): # The displayable image of a line.
    my_slope=line.left_of().slope(line.right_of())
    return "[" + line.left_of().image() + "," + \
           line.right_of().image() + \
           " slope=" + `my_slope` + \
           " y_int=" + `line.left_of().y_intercept(my_slope)` + "]";
  def intersection_lambdas (left_line, right_line,infinite):
    # Computes the intersection of two lines as parameters, of which the
first is used.
    # This is designed to work on all 2-D lines including vertical
lines.
    # t1 R1 + (1-t1)L1 = p = t2 R2 + (1-t2) L2 gives:
    # (R1-L1)t1 + (L2-R2) t2 + (L1-L2) = 0     gives simultaneously:
    # e t1 + c t2 + a = 0
    # f t1 + d t2 + b = 0
    assert(right_line.is_line());
    a=left_line.left_of().x_of()  - right_line.left_of().x_of();
    b=left_line.left_of().y_of()  - right_line.left_of().y_of();
    c=right_line.left_of().x_of() - right_line.right_of().x_of();
    d=right_line.left_of().y_of() - right_line.right_of().y_of();
    e=left_line.right_of().x_of() - left_line.left_of().x_of();
    f=left_line.right_of().y_of() - left_line.left_of().y_of();
    denominator=d*e-c*f;
    if close(denominator,0.0,line_threshold):
      return [infinite,infinite];
    else:  
      numerator=a*f-b*e;
      lambda2 = numerator/float(denominator);
    if   not close(e,0.0,line_threshold):  
      lambda1 = (-a-c*lambda2)/e; 
    elif not close(f,0.0,line_threshold):        
      lambda1=(-b-d*lambda2)/f
    else: 
      return [infinite,infinite];
    return [lambda1, lambda2];
  def is_line(line): # Type check.
    return 1;  
  def parameter_of(line,point): # Converts a point to parameter form.
    # Returns the parameter of a point on a line. Assumes that the point
is on the line.
    assert(point.is_point2D());
    numerator  =point-line.left_of();
    denominator=line.right_of()-line.left_of();
    if denominator.close(points2D(0.0,0.0),line_threshold):
      raise bad_line;
    if not close(denominator.x_of(),0.0,line_threshold): 
      return numerator.x_of()/denominator.x_of();
    else:
      return numerator.y_of()/denominator.y_of();
  def parameter2Point(line,llambda): # Converts a parameter to point
form.
    difference = line.right_of()-line.left_of();
    return line.left_of()+difference.times_scalar(llambda);
  def print_segment_intersection(line1,line2):
    t1,t2=line1.intersection_lambdas(line2,20000000.0)
    print "-----------------------------------------------"
    print "Line 1="+line1.image()
    print "Line 2="+line2.image()
    print "Intersection point for parameter +`t1` is " +
line1.parameter2Point(t1).image()
    print "The intersection point is on the segment: "
+["no","yes"][0.0<=t1<=1.0]
  
  
def linetest():
  P1=points2D(5,0)
  P2=points2D(0,5)
  L1=lines(P1,P2)
  P3=points2D(0,0)
  P4=points2D(1,1)
  L2=lines(P3,P4)
  L1.print_segment_intersection(L2)

linetest()