Function msDistanceSegmentToSegment

mapsearch.c:378–462

/* Users of this code must verify correctness for their application. */

double msDistanceSegmentToSegment(pointObj *pa, pointObj *pb, pointObj *pc, pointObj *pd) 
{
  vectorObj dP;
  vectorObj u, v, w;
  double a, b, c, d, e; 
  double D;
  double sc, sN, sD; /* N=numerator, D=demoninator */
  double tc, tN, tD;

  /* check for strictly parallel segments first */
  /* if(((pa->x == pb->x) && (pc->x == pd->x)) || (slope(pa,pb) == slope(pc,pd))) { // vertical (infinite slope) || otherwise parallel */
  /* D = msDistancePointToSegment(pa, pc, pd); */
  /* D = MS_MIN(D, msDistancePointToSegment(pb, pc, pd)); */
  /* D = MS_MIN(D, msDistancePointToSegment(pc, pa, pb)); */
  /* return(MS_MIN(D, msDistancePointToSegment(pd, pa, pb))); */
  /* } */

  u.x = pb->x - pa->x; /* u = pb - pa */
  u.y = pb->y - pa->y;
  v.x = pd->x - pc->x; /* v = pd - pc  */
  v.y = pd->y - pc->y; 
  w.x = pa->x - pc->x; /* w = pa - pc */
  w.y = pa->y - pc->y; 

  a = dot(u,u);
  b = dot(u,v);
  c = dot(v,v);
  d = dot(u,w);
  e = dot(v,w);

  D = a*c - b*b;
  sc = sN = sD = D;
  tc = tN = tD = D;
  
  /* compute the line parameters of the two closest points */
  if(D < SMALL_NUMBER) { /* lines are parallel or almost parallel */
    sN = 0.0;
    sD = 1.0;
    tN = e;
    tD = c;
  } else { /* get the closest points on the infinite lines */
    sN = b*e - c*d;
    tN = a*e - b*d;
    if(sN < 0) {
      sN = 0.0;
      tN = e;
      tD = c;
    } else if(sN > sD) {
      sN = sD;
      tN = e + b;
      tD = c;
    }
  }

  if(tN < 0) {
    tN = 0.0;
    if(-d < 0)
      sN = 0.0;