Skip to content

Instantly share code, notes, and snippets.

@aikar
Created August 13, 2010 14:43
Show Gist options
  • Save aikar/522991 to your computer and use it in GitHub Desktop.
Save aikar/522991 to your computer and use it in GitHub Desktop.
<?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;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment