Listing 1: The intersect function

int sgn(double x)
{
  if (abs(x) < mch_error) return 0;
  return (x<0.0) ? -1 : 1;
}
int triangulation::intersect(vector a, segment b)
{
  //  First check for zero or vertical vectors.
  if (eq_2d(a[0],a[1]) || eq_2d(b[0], b[1]))
    return true;

  // Now check that there are not points on top of other points.
  if ((!(a[0]==b[0]) && eq_2d(a[0],b[0])) ||
      (!(a[0]==b[1]) && eq_2d(a[0],b[1])) ||
      (!(a[1]==b[0]) && eq_2d(a[1],b[0])) ||
      (!(a[1]==b[1]) && eq_2d(a[1],b[1])))
    return true;

  //  Next check if the two vectors are the same.
  if ((eq_2d(a[0],b[0]) && eq_2d(a[1],b[1])) ||
      (eq_2d(a[0],b[1]) && eq_2d(a[1],b[0])))
    return true;

  int test1=sgn(dot_2d(perp_2d(a[1]-a[0]), b[0]-a[0]) *
    dot_2d(perp_2d(a[1]-a[0]), b[1]-a[0]));

  if (test1 > 0) return false;

  int test2=sgn(dot_2d(perp_2d(b[1]-b[0]), a[0]-b[0]) *
    dot_2d(perp_2d(b[1]-b[0]), a[1]-b[0]));

  if (test2 > 0) return false;
  if (test1<0 || test2<0)
    return true;

  //  if execution gets here, there are colinear points.
  //  Check to see if the two vectors are parallel.
  //  If they are not parallel then they share only one end point and
  //  are considered not to intersect.
  if (sgn(dot_2d(perp_2d(a[1]-a[0]), b[1]-b[0])) != 0) return false;

  //  if execution reaches here then the segments lie on the same line.
  //  All that is left is to determine if the two segments overlap at all.
  if (between(b[0],a[0],a[1]) || between(b[1],a[0],a[1]))
    return true;
  return between(a[0], b[0], b[1]);
}
//End of File