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?

Asked 7 Months ago Answers: 5 Viewed 49 times

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?

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:

**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; const double rad=180.0/M_PI; // 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 !!!**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:

**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`

.**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.**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 rad=180.0/M_PI;
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 pos0[3]={ 23.000000*deg, 48.000000*deg,2500.000000 }; // [rad,rad,m]
double pos1[3]={ 23.000000*deg, 35.000000*deg,2500.000000 }; // [rad,rad,m]
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;
}
//---------------------------------------------------------------------------
void WGS84_m2rad(double &da,double &db,double *abh)
{
// 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;
WGS84_m2rad(da,db,p0);
p0[0]+=vel0[0]*dt*da;
p0[1]+=vel0[1]*dt*db;
p0[2]+=vel0[2]*dt;
WGS84_m2rad(da,db,p1);
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

answered 4 Months ago

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

answered 4 Months ago

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

answered 4 Months ago

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

answered 4 Months ago

Only authorized users can answer the question. Please sign in first, or register a free account.

Not the answer you're looking for? Browse other questions tagged :

Taking

Eis the starting point of the ray,Lis the end point of the ray,Cis the center of sphere you're testing againstris the radius of that sphereCompute:

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 * dThis is a parametric equation:

P

_{x}= E_{x}+ td_{x}P

_{y}= E_{y}+ td_{y}into

(x - h)^{2}+ (y - k)^{2}= r^{2}(h,k) = center of circle.

to get:Expandx^{2}- 2xh + h^{2}+ y^{2}- 2yk + k^{2}- r^{2}= 0Plugx = e_{x}+ td_{x}y = e

_{y}+ td_{y}( e

_{x}+ td_{x})^{2}- 2( e_{x}+ td_{x})h + h^{2}+ ( e_{y}+ td_{y})^{2}- 2( e_{y}+ td_{y})k + k^{2}- r^{2}= 0Explodee_{x}^{2}+ 2e_{x}td_{x}+ t^{2}d_{x}^{2}- 2e_{x}h - 2td_{x}h + h^{2}+ e_{y}^{2}+ 2e_{y}td_{y}+ t^{2}d_{y}^{2}- 2e_{y}k - 2td_{y}k + k^{2}- r^{2}= 0Groupt^{2}( d_{x}^{2}+ d_{y}^{2}) + 2t( e_{x}d_{x}+ e_{y}d_{y}- d_{x}h - d_{y}k ) + e_{x}^{2}+ e_{y}^{2}- 2e_{x}h - 2e_{y}k + h^{2}+ k^{2}- r^{2}= 0Finally,t^{2}(d·d) + 2t(e·d-d·c) +e·e- 2(e·c) +c·c- r^{2}= 0Where

dis the vector d and · is the dot product.And then,t^{2}(d·d) + 2t(d· (e-c) ) + (e-c) · (e-c) - r^{2}= 0Lettingf=e-ct^{2}(d·d) + 2t(d·f) +f·f- r^{2}= 0So we get:

t

^{2}* (d·d) + 2t*(f·d) + (f·f- r^{2}) = 0So solving the quadratic equation: