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.

Image

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?

 Answers

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

So solving the quadratic equation:

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
 
CodeCaster
answered 7 Months ago
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;
    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 !!!

  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 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:

overview

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
 
cusejuice
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
 
turik
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
 
Shreejibawa
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
 
j_random_hacker
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 :  
Share