Find intersection of circle and line

Question

I am trying to find a line is intersecting a circle or not.

I write code but seems some issues with that code.

private Point2d[] IntersectionPoint(Point2d p1, Point2d p2, Point2d sc, double r)
    {
        Point2d[] sect = null;
        double a, b, c;
        double bb4ac;
        double mu1, mu2;
        Point2d dp;

        dp = p2 - p1;

        a = dp.X * dp.X + dp.Y * dp.Y;
        b = 2 * (dp.X * (p1.X - sc.X) + dp.Y * (p1.Y - sc.Y));
        c = sc.X * sc.X + sc.Y * sc.Y;
        c += p1.X * p1.X + p1.Y * p1.Y;
        c -= 2 * (sc.X * p1.X + sc.Y * p1.Y);
        c -= r * r;
        bb4ac = b * b - 4 * a * c;

        if (Math.Abs(a) < Double.Epsilon || bb4ac < 0)
        {
            return new Point2d[0];
        }

        mu1 = (-b + Math.Sqrt(bb4ac)) / (2 * a);
        mu2 = (-b - Math.Sqrt(bb4ac)) / (2 * a);

        // no intersection
        if ((mu1 < 0 || mu1 > 1) && (mu2 < 0 || mu2 > 1))
        {
            sect = new Point2d[0];
        }
        // one point on mu1
        else if (mu1 > 0 && mu1 < 1 && (mu2 < 0 || mu2 > 1))
        {
            sect = new Point2d[1];
            sect[0] = p1 + ((p2 - p1) * mu1);
        }
        // one point on mu2
        else if (mu2 > 0 && mu2 < 1 && (mu1 < 0 || mu1 > 1))
        {
            sect = new Point2d[1];
            sect[0] = p1 + ((p2 - p1) * mu2);
        }
        //  one or two points
        else if (mu1 > 0 && mu1 < 1 && mu2 > 0 && mu2 < 1)
        {
            //  tangential
            if (mu1 == mu2)
            {
                sect = new Point2d[1];
                sect[0] = p1 + ((p2 - p1) * mu1);
            }
            //  two points
            else
            {
                sect = new Point2d[2];
                sect[0] = p1 + ((p2 - p1) * mu1);
                sect[1] = p1 + ((p2 - p1) * mu2);
            }
        }
        else
        {
            //  should NEVER get here
            sect = new Point2d[0];
        }
        return sect;
    }

And calling this function like

Point ptOld = points[oldPoint];
Point ptNew = points[newPoint];
Point2d p1 = new Point2d((float)ptOld.latitude, (float)ptOld.longitude);
Point2d p2 = new Point2d((float)ptNew.latitude, (float)ptNew.longitude);
Point2d sc = new Point2d((float)loc.latitude, (float)loc.longitude);

It fails when i am trying with these co-ordinates

30,-30

80,-40

10

https://www.dropbox.com/s/38r9eylt2p4xfvw/graph.png

Answer 1

Edit: Vectors are now available in C#.

There is my code for intersecting:

public static Vector3? IntersectRayCircle(Vector3 rayStart, Vector3 rayPoint, Vector3 circlePosition, float circleRadiusSquared)
{
    if (rayStart == rayPoint || circleRadiusSquared <= 0)
    {
        return null;
    }

    Vector3 nearest = GetNearestPoint(circlePosition, rayStart, rayPoint, false, false);
    float distanceSquared = Vector3.DistanceSquared(nearest, circlePosition);

    if (distanceSquared > circleRadiusSquared)
    {
        return null;
    }

    Vector3 offset = Vector3.Normalize(rayPoint - rayStart) * (float)Math.Sqrt(circleRadiusSquared - distanceSquared);

    if (Vector3.DistanceSquared(circlePosition, rayStart) < circleRadiusSquared)
    {
        return nearest + offset;
    }
    else
    {
        return nearest - offset;
    }
}

---

public static Vector3 GetNearestPoint(Vector3 location, Vector3 segmentStart, Vector3 segmentEnd, bool trimStart, bool trimEnd)
{
    if (segmentStart == segmentEnd)
    {
        throw new ArgumentException("segmentStart cannot be equal to segmentEnd.");
    }

    Vector3 AP = location - segmentStart;
    Vector3 AB = segmentEnd - segmentStart;

    float magnitudeAB = AB.LengthSquared();
    float ABAPproduct = Vector3.Dot(AP, AB);
    float distance = ABAPproduct / magnitudeAB;

    return (distance < 0 && trimStart) ? segmentStart : (distance > 1 && trimEnd) ? segmentEnd : segmentStart + AB * distance;
}

Answer 2

You could do some linear algebra:

1. Express the line as an origin point P1 and a normalized direction vector N

2. Project the center C of the circle onto the line: PC = P1 + dot(C - P1, N) * N

3. Compute the distance squared dSquared between points C and PC.

4. If equal (with some small tolerance) to radiusSquared, PC is on the circle and is the single intersection point.

5. If greater than radiusSquared, no intersection.

6. Otherwise, the two intersections are given by

   1. offset = sqrt(radiusSquared - dSquared).
   2. Intersections = PC +/- offset * N.