0

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

Prince Chopra
  • 167
  • 11
  • http://paulbourke.net/geometry/circlesphere/raysphere.c // After reviewing the code it looks like that is actually where you base your code from. // I had based some of my code from there too, I had an issue that `bb4ac` sometimes got so close to zero (but still shouldn't be; tangent lines to a sphere) that it would fail the intersection test. // I think you need to compare `bb4ac` to epsilon – AnotherUser Aug 08 '14 at 12:36
  • I am not able to build this code. @AnotherUser – Prince Chopra Aug 08 '14 at 13:20
  • 3
    I also don't believe someone who writes code containing mathematical formulas, can't provide any information about what's going wrong besides "seems some issues" and "it fails". Either the code does not match your design (you can debug that) or your design is flawed and we should not be looking at code _yet_. – C.Evenhuis Aug 08 '14 at 13:49
  • You really are checking whether the radius is less than the distance of the middle from the line. [See here..](http://stackoverflow.com/questions/4438244/how-to-calculate-shortest-2d-distance-between-a-point-and-a-line-segment-in-all) – TaW Aug 08 '14 at 13:55
  • Yes i am following graph to check it. @TaW – Prince Chopra Aug 08 '14 at 14:05
  • I draw this graph. Please see this "https://www.dropbox.com/s/38r9eylt2p4xfvw/graph.png" @TaW – Prince Chopra Aug 08 '14 at 14:08
  • @PrinceChopra I did not link you to that for you to build it, not sure why you even tried to. Its there for you to review. I would almost bet that your code comes from that link after reviewing the variable names and basic structure. // You need to debug your code, you have not asked a question, you have asked us to debug your code. – AnotherUser Aug 08 '14 at 14:28
  • The link is dead. Also: You need to __clarify__: do you want to know __IF__ the line intersects the circle __or__ do you also need the intersection points?? – TaW Aug 08 '14 at 14:35

2 Answers2

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;
}
AgentFire
  • 8,944
  • 8
  • 43
  • 90
1

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.
dbc
  • 104,963
  • 20
  • 228
  • 340