intersect.php

<?php
function intersectPoint($line1, $angle1, $line2, $angle2)
{
  $lat1 = deg2rad($line1['lat']);
  $lon1 = deg2rad($line1['lng']);
  $lat2 = deg2rad($line2['lat']);
  $lon2 = deg2rad($line2['lng']);
  $brng13 = deg2rad($angle1);
  $brng23 = deg2rad($angle2);
  $dLat = $lat2-$lat1;
  $dLon = $lon2-$lon1;

  $dist12 = 2*asin( sqrt( sin($dLat/2)*sin($dLat/2) +
    cos($lat1)*cos($lat2)*sin($dLon/2)*sin($dLon/2) ) );
  if ($dist12 == 0) return null;

  // initial/final bearings between points
  $brngA = acos( ( sin($lat2) - sin($lat1)*cos($dist12) ) /
    ( sin($dist12)*cos($lat1) ) );
  if (!is_numeric($brngA)) $brngA = 0;  // protect against rounding
  $brngB = acos( ( sin($lat1) - sin($lat2)*cos($dist12) ) /
    ( sin($dist12)*cos($lat2) ) );

  if (sin(lon2-$lon1) > 0) {
    $brng12 = $brngA;
    $brng21 = 2*M_PI - $brngB;
  } else {
    $brng12 = 2*M_PI - $brngA;
    $brng21 = $brngB;
  }

  $alpha1 = ($brng13 - $brng12 + M_PI) % (2*M_PI) - M_PI;  // angle 2-1-3
  $alpha2 = ($brng21 - $brng23 + M_PI) % (2*M_PI) - M_PI;  // angle 1-2-3

  if (sin($alpha1)==0 && sin($alpha2)==0) return null;  // infinite intersections
  if (sin($alpha1)*sin($alpha2) < 0) return null;       // ambiguous intersection

  //alpha1 = abs(alpha1);
  //alpha2 = abs(alpha2);
  // ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?

  $alpha3 = acos( -cos($alpha1)*cos($alpha2) +
                       sin($alpha1)*sin($alpha2)*cos($dist12) );
  $dist13 = atan2( sin($dist12)*sin($alpha1)*sin($alpha2),
                       cos($alpha2)+cos($alpha1)*cos($alpha3) );
  $lat3 = asin( sin($lat1)*cos($dist13) +
                    cos($lat1)*sin($dist13)*cos($brng13) );
  $dLon13 = atan2( sin($brng13)*sin($dist13)*cos($lat1),
                       cos($dist13)-sin($lat1)*sin($lat3) );
  $lon3 = $lon1+$dLon13;
  $lon3 = ($lon3+M_PI) % (2*M_PI) - M_PI;  // normalise to -180..180º

  return array('lat' => rad2deg($lat3), 'lng' => rad2deg($lon3));
}
function intersectPoint($line1start, $line1end, $line2start, $line2end)
{
    $a1 = $line1start;
    $a2 = $line1end;
    $b1 = $line2start;
    $b2 = $line2end;

    $ua_t = ($b2['lat'] - $b1['lat']) * ($a1['lng'] - $b1['lng']) - ($b2['lng'] - $b1['lng']) * ($a1['lat'] - $b1['lat']);
    $ub_t = ($a2['lat'] - $a1['lat']) * ($a1['lng'] - $b1['lng']) - ($a2['lng'] - $a1['lng']) * ($a1['lat'] - $b1['lat']);
    $u_b  = ($b2['lng'] - $b1['lng']) * ($a2['lat'] - $a1['lat']) - ($b2['lat'] - $b1['lat']) * ($a2['lng'] - $a1['lng']);

    if ($u_b != 0 ) {
        $ua = $ua_t / $u_b;
        $ub = $ub_t / $u_b;

        if ( 0 <= $ua && $ua <= 1 && 0 <= $ub && $ub <= 1 ) {
            return array(
                'lat' => $a1['lat'] + $ua * ($a2['lat'] - $a1['lat']),
                'lng' => $a1['lng'] + $ua * ($a2['lng'] - $a1['lng'])
                
            );
        }
    }

    return null;
};
function intersectPoint($line1start, $line1end, $line2start, $line2end)
//($p0_x, $p0_y, $p1_x, $p1_y, $p2_x, $p2_y, $p3_x, $p3_y)
{
    $p0_x = $line1start['lat'];
    $p0_y = $line1start['lng'];
    $p1_x = $line1end['lat'];
    $p1_y = $line1end['lng'];
    $p3_x = $line2start['lat'];
    $p3_y = $line2start['lng'];
    $p4_x = $line1end['lat'];
    $p4_y = $line1end['lng'];
    
    $s1_x = $p1_x - $p0_x;     $s1_y = $p1_y - $p0_y;
    $s2_x = $p3_x - $p2_x;     $s2_y = $p3_y - $p2_y;


    $s = (-$s1_y * ($p0_x - $p2_x) + $s1_x * ($p0_y - $p2_y)) / (-$s2_x * $s1_y + $s1_x * $s2_y);
    $t = ( $s2_x * ($p0_y - $p2_y) - $s2_y * ($p0_x - $p2_x)) / (-$s2_x * $s1_y + $s1_x * $s2_y);

    if ($s >= 0 && $$s <= 1 && $t >= 0 && $t <= 1)
    {
        // Colli$sion detected
        return array(
            'lat' => $p0_x + ($t * $s1_x),
            'lng' => $p0_y + ($t * $s1_y)
        );
        return 1;
    }

    return null; // No collision
}
class GEO
{
   const
    PI                = M_PI,
    TWOPI             = 6.28318530718,
    DE2RA             = 0.01745329252,
    RA2DE             = 57.2957795129,
    ERAD              = 6378.135,
    ERADM             = 6378135.0,
    AVG_ERAD          = 6371.0,
    FLATTENING        = 0.00335281066,
                                  // Earth flattening (WGS '84)
    EPS               = 0.000000000005,
    KM2MI             = 0.621371,
    GEOSTATIONARY_ALT = 35786.0;    // km
}

function intersectPoint($line1start, $line1end, $line2start, $line2end)
{
    $lat1A = $line1start['lat'];
    $lon1A = $line1start['lng'];
    $lat1B = $line1end['lat'];
    $lon1B = $line1end['lng'];

    $lat2A = $line2start['lat'];
    $lon2A = $line2start['lng'];
    $lat2B = $line2end['lat'];
    $lon2B = $line2end['lng'];

   $isInside = false;

   $v2 = $v1 = array();

   

   $d1 = distance($lat1A, $lon1A, $lat1B, $lon1B);
   $d2 = distance($lat2A, $lon2A, $lat2B, $lon2B);

   //
   // for path 1, setting up my 2 vectors, $v1 is vector
   // from center of the Earth to point A, $v2 is vector
   // from center of the Earth to point B.
   //
   $v1[0] = cos(deg2rad($lat1A)) * cos(deg2rad($lon1A));
   $v1[1] = sin(deg2rad($lat1A));
   $v1[2] = cos(deg2rad($lat1A)) * sin(deg2rad($lon1A));

   $v2[0] = cos(deg2rad($lat1B)) * cos(deg2rad($lon1B));
   $v2[1] = sin(deg2rad($lat1B));
   $v2[2] = cos(deg2rad($lat1B)) * sin(deg2rad($lon1B));

   //
   // $n1 is the normal to the plane formed by the three points:
   // center of the Earth, point 1A, and point 1B
   //
   $n1[0] = ($v1[1]*$v2[2]) - ($v1[2]*$v2[1]);
   $n1[1] = ($v1[2]*$v2[0]) - ($v1[0]*$v2[2]);
   $n1[2] = ($v1[0]*$v2[1]) - ($v1[1]*$v2[0]);

   //
   // for path 2, setting up my 2 vectors, $v1 is vector
   // from center of the Earth to point A, $v2 is vector
   // from center of the Earth to point B.
   //
   $v1[0] = cos(deg2rad($lat2A)) * cos(deg2rad($lon2A));
   $v1[1] = sin(deg2rad($lat2A));
   $v1[2] = cos(deg2rad($lat2A)) * sin(deg2rad($lon2A));

   $v2[0] = cos(deg2rad($lat2B)) * cos(deg2rad($lon2B));
   $v2[1] = sin(deg2rad($lat2B));
   $v2[2] = cos(deg2rad($lat2B)) * sin(deg2rad($lon2B));

   //
   // $n1 is the normal to the plane formed by the three points:
   // center of the Earth, point 2A, and point 2B
   //
   $n2[0] = ($v1[1]*$v2[2]) - ($v1[2]*$v2[1]);
   $n2[1] = ($v1[2]*$v2[0]) - ($v1[0]*$v2[2]);
   $n2[2] = ($v1[0]*$v2[1]) - ($v1[1]*$v2[0]);

   //
   // $v3 is perpendicular to both normal 1 and normal 2, so
   // it lies in both planes, so it must be the line of
   // intersection of the planes. The question is: does it
   // go towards the correct intersection point or away from
   // it.
   //
   $v3a[0] = ($n1[1]*$n2[2]) - ($n1[2]*$n2[1]);
   $v3a[1] = ($n1[2]*$n2[0]) - ($n1[0]*$n2[2]);
   $v3a[2] = ($n1[0]*$n2[1]) - ($n1[1]*$n2[0]);

   //
   // want to make $v3a a unit vector, so first have to get
   // magnitude, then each component by magnitude
   //
   $m = sqrt($v3a[0]*$v3a[0] + $v3a[1]*$v3a[1] + $v3a[2]*$v3a[2]);
   $v3a[0] /= $m;
   $v3a[1] /= $m;
   $v3a[2] /= $m;

   //
   // calcu$lating intersection points 3A & 3B. A & B are
   // arbitrary designations right now, $later we make A
   // the one close to, or within, the paths.
   //
   $lat3A = asin($v3a[1]);
   if (($lat3A > GEO::EPS) || (-$lat3A > GEO::EPS))
      $lon3A = asin($v3a[2]/cos($lat3A));
   else
      $lon3A = 0.0;

   $v3b[0] = ($n2[1]*$n1[2]) - ($n2[2]*$n1[1]);
   $v3b[1] = ($n2[2]*$n1[0]) - ($n2[0]*$n1[2]);
   $v3b[2] = ($n2[0]*$n1[1]) - ($n2[1]*$n1[0]);

   $m = sqrt($v3b[0]*$v3b[0] + $v3b[1]*$v3b[1] + $v3b[2]*$v3b[2]);
   $v3b[0] /= $m;
   $v3b[1] /= $m;
   $v3b[2] /= $m;

   $lat3B = asin($v3b[1]);
   if (($lat3B > GEO::EPS) || (-$lat3B > GEO::EPS))
      $lon3B = asin($v3b[2]/cos($lat3B));
   else
      $lon3B = 0.0;

   //
   // get the distances from the path endpoints to the two
   // intersection points. these values will be used to determine
   // which intersection point lies on the paths, or which one
   // is closer.
   //
   $d1a3a = distance($lat1A, $lon1A, $lat3A, $lon3A);
   $d1b3a = distance($lat1B, $lon1B, $lat3A, $lon3A);
   $d2a3a = distance($lat2A, $lon2A, $lat3A, $lon3A);
   $d2b3a = distance($lat2B, $lon2B, $lat3A, $lon3A);

   $d1a3b = distance($lat1A, $lon1A, $lat3B, $lon3B);
   $d1b3b = distance($lat1B, $lon1B, $lat3B, $lon3B);
   $d2a3b = distance($lat2A, $lon2A, $lat3B, $lon3B);
   $d2b3b = distance($lat2B, $lon2B, $lat3B, $lon3B);

   if (($d1a3a < $d1) && ($d1b3a < $d1) && ($d2a3a < $d2)
                    && ($d2b3a < $d2))
   {
      $isInside = true;
   }
   else if (($d1a3b < $d1) && ($d1b3b < $d1) && ($d2a3b < $d2)
                         && ($d2b3b < $d2))
   {
      //
      // 3b is inside the two paths, so swap 3a & 3b
      //
      $isInside = true;

      $m = $lat3A;
      $lat3A = $lat3B;
      $lat3B = $m;
      $m = $lon3A;
      $lon3A = $lon3B;
      $lon3B = $m;
   }
   else
   {
      //
      // figure out which one is closer to the path
      //
      $d1 = $d1a3a + $d1b3a + $d2a3a + $d2b3a;
      $d2 = $d1a3b + $d1b3b + $d2a3b + $d2b3b;

      if ($d1 > $d2)
      {
         //
         // Okay, we are here because 3b {$lat3B,$lon3B} is closer to
         // the paths, so we need to swap 3a & 3b. the other case is
         // already the way 3a & 3b are organized, so no need to swap
         //
         $m = $lat3A;
         $lat3A = $lat3B;
         $lat3B = $m;
         $m = $lon3A;
         $lon3A = $lon3B;
         $lon3B = $m;
      }
   }

   //
   // convert the intersection points to degrees
   //
   $lat3A *= GEO::RA2DE;
   $lon3A *= GEO::RA2DE;
   $lat3B *= GEO::RA2DE;
   $lon3B *= GEO::RA2DE;

   if($lat3A && $lon3A)
   {
       return array('lat' => $lat3A, 'lng' => $lon3A, 'e' => array($lat3B, $lon3B), 'i' => ($isInside ? 'true' : 'false'));
   }else{
       return null;
   }
   
}