# Circle line-segment collision detection algorithm?

I have a line from A to B and a circle positioned at C with the radius R.

What is a good algorithm to use to check whether the line intersects the circle? And at what coordinate along the circles edge it occurred?

60

Taking

1. E is the starting point of the ray,
2. L is the end point of the ray,
3. C is the center of sphere you're testing against
4. r is the radius of that sphere

Compute:
d = L - E ( Direction vector of ray, from start to end )
f = E - C ( Vector from center sphere to ray start )

Then the intersection is found by..
Plugging:
P = E + t * d
This is a parametric equation:
Px = Ex + tdx
Py = Ey + tdy
into
(x - h)2 + (y - k)2 = r2
(h,k) = center of circle.

Note: We've simplified the problem to 2D here, the solution we get applies also in 3D

to get:

1. Expand x2 - 2xh + h2 + y2 - 2yk + k2 - r2 = 0
2. Plug x = ex + tdx
y = ey + tdy
( ex + tdx )2 - 2( ex + tdx )h + h2 + ( ey + tdy )2 - 2( ey + tdy )k + k2 - r2 = 0
3. Explode ex2 + 2extdx + t2dx2 - 2exh - 2tdxh + h2 + ey2 + 2eytdy + t2dy2 - 2eyk - 2tdyk + k2 - r2 = 0
4. Group t2( dx2 + dy2 ) + 2t( exdx + eydy - dxh - dyk ) + ex2 + ey2 - 2exh - 2eyk + h2 + k2 - r2 = 0
5. Finally, t2( d · d ) + 2t( e · d - d · c ) + e · e - 2( e · c ) + c · c - r2 = 0
Where d is the vector d and · is the dot product.
6. And then, t2( d · d ) + 2t( d · ( e - c ) ) + ( e - c ) · ( e - c ) - r2 = 0
7. Letting f = e - c t2( d · d ) + 2t( d · f ) + f · f - r2 = 0

So we get:
t2 * (d · d) + 2t*( f · d ) + ( f · f - r2 ) = 0

``````float a = d.Dot( d ) ;
float b = 2*f.Dot( d ) ;
float c = f.Dot( f ) - r*r ;

float discriminant = b*b-4*a*c;
if( discriminant < 0 )
{
// no intersection
}
else
{
// ray didn't totally miss sphere,
// so there is a solution to
// the equation.

discriminant = sqrt( discriminant );

// either solution may be on or off the ray so need to test both
// t1 is always the smaller value, because BOTH discriminant and
// a are nonnegative.
float t1 = (-b - discriminant)/(2*a);
float t2 = (-b + discriminant)/(2*a);

// 3x HIT cases:
//          -o->             --|-->  |            |  --|->
// Impale(t1 hit,t2 hit), Poke(t1 hit,t2>1), ExitWound(t1<0, t2 hit),

// 3x MISS cases:
//       ->  o                     o ->              | -> |
// FallShort (t1>1,t2>1), Past (t1<0,t2<0), CompletelyInside(t1<0, t2>1)

if( t1 >= 0 && t1 <= 1 )
{
// t1 is the intersection, and it's closer than t2
// (since t1 uses -b - discriminant)
// Impale, Poke
return true ;
}

// here t1 didn't intersect so we are either started
// inside the sphere or completely past it
if( t2 >= 0 && t2 <= 1 )
{
// ExitWound
return true ;
}

// no intn: FallShort, Past, CompletelyInside
return false ;
}
``````
Tuesday, June 1, 2021

31

[Edit5] Complete reedit in case you need old sources see the revision history

As Nico Schertler pointed out checking for line to line intersection is insanity as the probability of intersecting 2 trajectories at same position and time is almost none (even when including round-off precision overlaps). Instead you should find place on each trajectory that is close enough (to collide) and both objects are there at relatively same time. Another problem is that your trajectories are not linear at all. Yes they can appear linear for shor times in booth WGS84 and Cartesian but with increasing time the trajectory bends around Earth. Also your input values units for speed are making this a bit harder so let me recapitulate normalized values I will be dealing with from now:

1. Input

consists of two objects. For each is known its actual position (in WGS84 `[rad]`) and actual speeds `[m/s]` but not in Cartesian space but WGS84 local axises instead. For example something like this:

``````const double kmh=1.0/3.6;
const double deg=M_PI/180.0;
//                      lon            lat      alt
double pos0[3]={  23.000000*deg, 48.000000*deg,2500.000000 };
double pos1[3]={  23.000000*deg, 35.000000*deg,2500.000000 };
double vel0[3]={ 100.000000*kmh,-20.000000*kmh,   0.000000*kmh };
double vel1[3]={ 100.000000*kmh, 20.000000*kmh,   0.000000*kmh };
``````

Beware mine coordinates are in `Long,Lat,Alt` order/convention !!!

2. output

You need to compute the time in which the two objects "collide" Additional constrains to solution can be added latter on. As mentioned before we are not searching for intersection but "closest" approach instead that suffice collision conditions (like distance is smaller then some threshold).

After some taught and testing I decided to use iterative approach in WGS84 space. That brings up some problems like how to convert speed in `[m/s]` in WGS84 space to `[rad/s]` in WGS84 space. This ratio is changing with object altitude and latitude. In reality we need to compute angle change in `long` and `lat` axises that are "precisely" equal to `1m` traveled distance and then multiply the velocities by it. This can be approximated by arc-length equation:

``````l = dang.R
``````

Where `R` is actual radius of angular movement, `ang` is the angle change and `l` is traveled distance so when `l=1.0` then:

``````dang = 1.0/R
``````

If we have Cartesian position `x,y,z` (`z` is Earth rotation axis) then:

``````Rlon = sqrt (x*x + y*y)
Rlat = sqrt (x*x + y*y + z*z)
``````

Now we can iterate positions with time which can be used to approximate closest approach time. We need to limit the max time step however so we do not miss to much of the Earth curvature. This limit is dependent on used speeds and target precision. So here the algorithm to find the approach:

1. init

set initial time step to the upper limit like `dt=1000.0` and compute actual positions of booth objects in Cartesian space. From that compute their distance `d1`.

2. iteration

set `d0=d1` then compute actual speeds in WGS84 for actual positions and add `speed*dt` to each objects actual `WGS84` position. Now just compute actual positions in Cartesian space and compute their distance `d1`

if `d0>d1` then it menas we are closing to the closest approach so goto #2 again.
In case `d0==d1` the trajectories are parallel so return approach time `t=0.0`
In case `d0<d1` we already crossed the closest approach so set `dt = -0.1*dt` and if `dt>=desired_accuracy` goto #2 otherwise stop.

3. recover best `t`

After the iteration in #2 we should recover the best time back so return `t+10.0*dt;`

Now we have closest approach time `t`. Beware it can be negative (if the objects are going away from each other). Now you can add your constrains like

``````if (d0<_max_d)
if ((t>=0.0)&&(t<=_max_T))
return collision ...
``````

Here C++ source for this:

``````//---------------------------------------------------------------------------
#include <math.h>
//---------------------------------------------------------------------------
const double kmh=1.0/3.6;
const double deg=M_PI/180.0;
const double  _earth_a=6378137.00000;   // [m] WGS84 equator radius
const double  _earth_b=6356752.31414;   // [m] WGS84 epolar radius
const double  _earth_e=8.1819190842622e-2; //  WGS84 eccentricity
const double  _earth_ee=_earth_e*_earth_e;
//--------------------------------------------------------------------------
const double _max_d=2500.0;         // [m] collision gap
const double _max_T=3600000.0;      // [s] max collision time
const double _max_dt=1000.0;        // [s] max iteration time step (for preserving accuracy)
//--------------------------------------------------------------------------
//                      lon            lat      alt
double vel0[3]={ 100.000000*kmh,-20.000000*kmh,   0.000000*kmh }; // [m/s,m/s,m/s]
double vel1[3]={ 100.000000*kmh,+20.000000*kmh,   0.000000*kmh }; // [m/s,m/s,m/s]
//---------------------------------------------------------------------------
double divide(double x,double y)
{
if ((y>=-1e-30)&&(y<=+1e-30)) return 0.0;
return x/y;
}
void  vector_copy(double *c,double *a)         { for(int i=0;i<3;i++) c[i]=a[i];       }
double vector_len(double *a) { return sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])); }
void  vector_len(double *c,double *a,double l)
{
l=divide(l,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
c[0]=a[0]*l;
c[1]=a[1]*l;
c[2]=a[2]*l;
}
void  vector_sub(double *c,double *a,double *b) { for(int i=0;i<3;i++) c[i]=a[i]-b[i]; }
//---------------------------------------------------------------------------
void WGS84toXYZ(double *xyz,double *abh)
{
double  a,b,h,l,c,s;
a=abh[0];
b=abh[1];
h=abh[2];
c=cos(b);
s=sin(b);
// WGS84 from eccentricity
l=_earth_a/sqrt(1.0-(_earth_ee*s*s));
xyz[0]=(l+h)*c*cos(a);
xyz[1]=(l+h)*c*sin(a);
xyz[2]=(((1.0-_earth_ee)*l)+h)*s;
}
//---------------------------------------------------------------------------
{
// WGS84 from eccentricity
double p[3],rr;
WGS84toXYZ(p,abh);
rr=(p[0]*p[0])+(p[1]*p[1]);
da=divide(1.0,sqrt(rr));
rr+=p[2]*p[2];
db=divide(1.0,sqrt(rr));
}
//---------------------------------------------------------------------------
double collision(double *pos0,double *vel0,double *pos1,double *vel1)
{
int e,i,n;
double p0[3],p1[3],q0[3],q1[3],da,db,dt,t,d0,d1,x,y,z;
vector_copy(p0,pos0);
vector_copy(p1,pos1);
// find closest d1[m] approach in time t[sec]
dt=_max_dt; // [sec] initial time step (accuracy = dt/10^(n-1)
n=6;        // acuracy loops
for (t=0.0,i=0;i<n;i++)
for (e=0;;e=1)
{
d0=d1;
// compute xyz distance
WGS84toXYZ(q0,p0);
WGS84toXYZ(q1,p1);
vector_sub(q0,q0,q1);
d1=vector_len(q0);
// nearest approach crossed?
if (e)
{
if (d0<d1){ dt*=-0.1; break; }                  // crossing trajectories
if (fabs(d0-d1)<=1e-10) { i=n; t=0.0; break; }  // parallel trajectories
}
// apply time step
t+=dt;
p0[0]+=vel0[0]*dt*da;
p0[1]+=vel0[1]*dt*db;
p0[2]+=vel0[2]*dt;
p1[0]+=vel1[0]*dt*da;
p1[1]+=vel1[1]*dt*db;
p1[2]+=vel1[2]*dt;
}
t+=10.0*dt; // recover original t
//  if ((d0<_max_d)&&(t>=0.0)&&(t<=_max_T)) return collision; else return no_collision;
return t;
}
//---------------------------------------------------------------------------
``````

Here an overview of example:

Red is object0 and Green is object1. The White squares represents position at computed collision at time `t_coll [s]` with distance `d_coll [m]`. Yellow squares are positions at user defined time `t_anim [s]` with distance `d_anim [m]` which is controlled by User for debugging purposes. As you can see this approach works also for times like 36 hours ...

Hope I did not forget to copy something (if yes comment me and I will add it)

Saturday, July 31, 2021

30

I could be wrong about this, but there are a few spots in the code that seem very suspicious. To begin, consider this line:

``````// calculate plane
float d = Dot(normal, coord);
``````

Here, your value `d` corresponds to the dot product between the plane normal (a vector) and a point in space (a point on the plane). This seems wrong. In particular, if you have any plane passing through the origin and use the origin as the coordinate point, you will end up computing

``````d = Dot(normal, (0, 0, 0)) = 0
``````

And immediately returning false. I'm not sure what you intended to do here, but I'm pretty sure that this isn't what you meant.

Another spot in the code that seems suspicious is this line:

``````// Compute the t value for the directed line ray intersecting the plane
float t = (d - Dot(normal, rayOrigin)) / Dot(normal, ray);
``````

Note that you're computing the dot product between the plane's normal vector (a vector) and the ray's origin point (a point in space). This seems weird because it means that depending on where the ray originates in space, the scaling factor you use for the ray changes. I would suggest looking at this code one more time to see if this is really what you meant.

Hope this helps!

Monday, August 9, 2021

53

A. Check if it is intersecting the whole cirlce.
B. Check if it is intersecting either of the straight segment lines.
C. If not, check if the angle between the circle centres lies in the angular range of the segment (dot product is good for this).

Intersection requires `A && (B || C)`

Tuesday, August 10, 2021

62

I don't see how you can represent a rectangle with just a "focus point". You'll need either the two corner points or one corner point with a width/height/rotation set of data.

However, once you have a rectangle, I would simply break it down into four lines and do the intercept checks between each of those four lines and the line you want to check against.

Doing a search on SO for "line intersection" turns up numerous questions, including this one, which seems promising. In fact, searching for "line rectangle intersection" gives you this one, which seems to be exactly what you're after.

Tuesday, August 24, 2021