我有一条从a到B的直线和一个半径为R的圆。

用什么算法来检查直线是否与圆相交?它在圆边的哪个坐标上?


当前回答

如果你找到了圆心(因为它是3D的,我想你是指球体而不是圆)和直线之间的距离,然后检查这个距离是否小于可以做到这一点的半径。

碰撞点显然是直线和球面之间最近的点(当你计算球面和直线之间的距离时,会计算出这个点)

点与线之间的距离: http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html

其他回答

好吧,我不会给你代码,但既然你已经标记了这个算法,我认为这对你来说无关紧要。 首先,你要得到一个垂直于这条直线的向量。

y = ax + c是一个未知变量c是未知变量 为了解决这个问题,计算直线经过圆心时的值。

也就是说, 将圆心的位置代入直线方程,解出c。 然后计算原直线与其法线的交点。

这样就能得到直线上离圆最近的点。 计算该点到圆中心之间的距离(使用矢量的大小)。 如果这个小于圆的半径,看,我们有一个交点!

我发现这个解决方案似乎比其他一些解决方案更容易遵循。

采取:

p1 and p2 as the points for the line, and
c as the center point for the circle and r for the radius

我可以用斜截式来解直线方程。但是,我不想处理以c为点的复杂方程,所以我只是平移了坐标系使圆在(0,0)处

p3 = p1 - c
p4 = p2 - c

顺便说一下,当我相互减分的时候,我是在减去x再减去y,然后把它们放到一个新的点里,以防有人不知道。

不管怎样,我现在解出p3和p4的直线方程

m = (p4_y - p3_y) / (p4_x - p3) (the underscore is an attempt at subscript)
y = mx + b
y - mx = b (just put in a point for x and y, and insert the m we found)

好的。现在我需要让这两个方程相等。首先我需要解圆的x方程

x^2 + y^2 = r^2
y^2 = r^2 - x^2
y = sqrt(r^2 - x^2)

然后我让它们相等:

mx + b = sqrt(r^2 - x^2)

求二次方程(0 = ax^2 + bx + c)

(mx + b)^2 = r^2 - x^2
(mx)^2 + 2mbx + b^2 = r^2 - x^2
0 = m^2 * x^2 + x^2 + 2mbx + b^2 - r^2
0 = (m^2 + 1) * x^2 + 2mbx + b^2 - r^2

现在我有了a b c。

a = m^2 + 1
b = 2mb
c = b^2 - r^2

我把这个代入二次公式

(-b ± sqrt(b^2 - 4ac)) / 2a

用值代入,然后尽可能简化:

(-2mb ± sqrt(b^2 - 4ac)) / 2a
(-2mb ± sqrt((-2mb)^2 - 4(m^2 + 1)(b^2 - r^2))) / 2(m^2 + 1)
(-2mb ± sqrt(4m^2 * b^2 - 4(m^2 * b^2 - m^2 * r^2 + b^2 - r^2))) / 2m^2 + 2
(-2mb ± sqrt(4 * (m^2 * b^2 - (m^2 * b^2 - m^2 * r^2 + b^2 - r^2))))/ 2m^2 + 2
(-2mb ± sqrt(4 * (m^2 * b^2 - m^2 * b^2 + m^2 * r^2 - b^2 + r^2)))/ 2m^2 + 2
(-2mb ± sqrt(4 * (m^2 * r^2 - b^2 + r^2)))/ 2m^2 + 2
(-2mb ± sqrt(4) * sqrt(m^2 * r^2 - b^2 + r^2))/ 2m^2 + 2
(-2mb ± 2 * sqrt(m^2 * r^2 - b^2 + r^2))/ 2m^2 + 2
(-2mb ± 2 * sqrt(m^2 * r^2 + r^2 - b^2))/ 2m^2 + 2
(-2mb ± 2 * sqrt(r^2 * (m^2 + 1) - b^2))/ 2m^2 + 2

这几乎是化简的极限了。最后,分离出带有±的方程:

(-2mb + 2 * sqrt(r^2 * (m^2 + 1) - b^2))/ 2m^2 + 2 or     
(-2mb - 2 * sqrt(r^2 * (m^2 + 1) - b^2))/ 2m^2 + 2 

然后简单地将这两个方程的结果代入mx + b中的x。为了清晰起见,我写了一些JavaScript代码来演示如何使用这个:

function interceptOnCircle(p1,p2,c,r){
    //p1 is the first line point
    //p2 is the second line point
    //c is the circle's center
    //r is the circle's radius

    var p3 = {x:p1.x - c.x, y:p1.y - c.y} //shifted line points
    var p4 = {x:p2.x - c.x, y:p2.y - c.y}

    var m = (p4.y - p3.y) / (p4.x - p3.x); //slope of the line
    var b = p3.y - m * p3.x; //y-intercept of line

    var underRadical = Math.pow((Math.pow(r,2)*(Math.pow(m,2)+1)),2)-Math.pow(b,2)); //the value under the square root sign 

    if (underRadical < 0){
    //line completely missed
        return false;
    } else {
        var t1 = (-2*m*b+2*Math.sqrt(underRadical))/(2 * Math.pow(m,2) + 2); //one of the intercept x's
        var t2 = (-2*m*b-2*Math.sqrt(underRadical))/(2 * Math.pow(m,2) + 2); //other intercept's x
        var i1 = {x:t1,y:m*t1+b} //intercept point 1
        var i2 = {x:t2,y:m*t2+b} //intercept point 2
        return [i1,i2];
    }
}

我希望这能有所帮助!

附注:如果任何人发现任何错误或有任何建议,请评论。我是新手,欢迎大家的帮助/建议。

以下是我在TypeScript中的解决方案,遵循@Mizipzor建议的想法(使用投影):

/**
 * Determines whether a line segment defined by a start and end point intersects with a sphere defined by a center point and a radius
 * @param a the start point of the line segment
 * @param b the end point of the line segment
 * @param c the center point of the sphere
 * @param r the radius of the sphere
 */
export function lineSphereIntersects(
  a: IPoint,
  b: IPoint,
  c: IPoint,
  r: number
): boolean {
  // find the three sides of the triangle formed by the three points
  const ab: number = distance(a, b);
  const ac: number = distance(a, c);
  const bc: number = distance(b, c);

  // check to see if either ends of the line segment are inside of the sphere
  if (ac < r || bc < r) {
    return true;
  }

  // find the angle between the line segment and the center of the sphere
  const numerator: number = Math.pow(ac, 2) + Math.pow(ab, 2) - Math.pow(bc, 2);
  const denominator: number = 2 * ac * ab;
  const cab: number = Math.acos(numerator / denominator);

  // find the distance from the center of the sphere and the line segment
  const cd: number = Math.sin(cab) * ac;

  // if the radius is at least as long as the distance between the center and the line
  if (r >= cd) {
    // find the distance between the line start and the point on the line closest to
    // the center of the sphere
    const ad: number = Math.cos(cab) * ac;
    // intersection occurs when the point on the line closest to the sphere center is
    // no further away than the end of the line
    return ad <= ab;
  }
  return false;
}

export function distance(a: IPoint, b: IPoint): number {
  return Math.sqrt(
    Math.pow(b.z - a.z, 2) + Math.pow(b.y - a.y, 2) + Math.pow(b.x - a.x, 2)
  );
}

export interface IPoint {
  x: number;
  y: number;
  z: number;
}

这里是一个用golang写的解决方案。这个方法和这里发布的其他一些答案类似,但不完全相同。它易于实现,并已经过测试。以下是步骤:

Translate coordinates so that the circle is at the origin. Express the line segment as parametrized functions of t for both the x and y coordinates. If t is 0, the function's values are one end point of the segment, and if t is 1, the function's values are the other end point. Solve, if possible, the quadratic equation resulting from constraining values of t that produce x, y coordinates with distances from the origin equal to the circle's radius. Throw out solutions where t is < 0 or > 1 ( <= 0 or >= 1 for an open segment). Those points are not contained in the segment. Translate back to original coordinates.

这里导出了二次曲线的A、B和C的值,其中(n-et)和(m-dt)分别是直线x坐标和y坐标的方程。R是圆的半径。

(n-et)(n-et) + (m-dt)(m-dt) = rr
nn - 2etn + etet + mm - 2mdt + dtdt = rr
(ee+dd)tt - 2(en + dm)t + nn + mm - rr = 0

因此A = ee+dd, B = - 2(en + dm), C = nn + mm - rr。

下面是函数的golang代码:

package geom

import (
    "math"
)

// SegmentCircleIntersection return points of intersection between a circle and
// a line segment. The Boolean intersects returns true if one or
// more solutions exist. If only one solution exists, 
// x1 == x2 and y1 == y2.
// s1x and s1y are coordinates for one end point of the segment, and
// s2x and s2y are coordinates for the other end of the segment.
// cx and cy are the coordinates of the center of the circle and
// r is the radius of the circle.
func SegmentCircleIntersection(s1x, s1y, s2x, s2y, cx, cy, r float64) (x1, y1, x2, y2 float64, intersects bool) {
    // (n-et) and (m-dt) are expressions for the x and y coordinates
    // of a parameterized line in coordinates whose origin is the
    // center of the circle.
    // When t = 0, (n-et) == s1x - cx and (m-dt) == s1y - cy
    // When t = 1, (n-et) == s2x - cx and (m-dt) == s2y - cy.
    n := s2x - cx
    m := s2y - cy

    e := s2x - s1x
    d := s2y - s1y

    // lineFunc checks if the  t parameter is in the segment and if so
    // calculates the line point in the unshifted coordinates (adds back
    // cx and cy.
    lineFunc := func(t float64) (x, y float64, inBounds bool) {
        inBounds = t >= 0 && t <= 1 // Check bounds on closed segment
        // To check bounds for an open segment use t > 0 && t < 1
        if inBounds { // Calc coords for point in segment
            x = n - e*t + cx
            y = m - d*t + cy
        }
        return
    }

    // Since we want the points on the line distance r from the origin,
    // (n-et)(n-et) + (m-dt)(m-dt) = rr.
    // Expanding and collecting terms yeilds the following quadratic equation:
    A, B, C := e*e+d*d, -2*(e*n+m*d), n*n+m*m-r*r

    D := B*B - 4*A*C // discriminant of quadratic
    if D < 0 {
        return // No solution
    }
    D = math.Sqrt(D)

    var p1In, p2In bool
    x1, y1, p1In = lineFunc((-B + D) / (2 * A)) // First root
    if D == 0.0 {
        intersects = p1In
        x2, y2 = x1, y1
        return // Only possible solution, quadratic has one root.
    }

    x2, y2, p2In = lineFunc((-B - D) / (2 * A)) // Second root

    intersects = p1In || p2In
    if p1In == false { // Only x2, y2 may be valid solutions
        x1, y1 = x2, y2
    } else if p2In == false { // Only x1, y1 are valid solutions
        x2, y2 = x1, y1
    }
    return
}

我用这个函数进行了测试,确认解点在线段内和圆上。它创建了一个测试段,并围绕给定的圆进行扫描:

package geom_test

import (
    "testing"

    . "**put your package path here**"
)

func CheckEpsilon(t *testing.T, v, epsilon float64, message string) {
    if v > epsilon || v < -epsilon {
        t.Error(message, v, epsilon)
        t.FailNow()
    }
}

func TestSegmentCircleIntersection(t *testing.T) {
    epsilon := 1e-10      // Something smallish
    x1, y1 := 5.0, 2.0    // segment end point 1
    x2, y2 := 50.0, 30.0  // segment end point 2
    cx, cy := 100.0, 90.0 // center of circle
    r := 80.0

    segx, segy := x2-x1, y2-y1

    testCntr, solutionCntr := 0, 0

    for i := -100; i < 100; i++ {
        for j := -100; j < 100; j++ {
            testCntr++
            s1x, s2x := x1+float64(i), x2+float64(i)
            s1y, s2y := y1+float64(j), y2+float64(j)

            sc1x, sc1y := s1x-cx, s1y-cy
            seg1Inside := sc1x*sc1x+sc1y*sc1y < r*r
            sc2x, sc2y := s2x-cx, s2y-cy
            seg2Inside := sc2x*sc2x+sc2y*sc2y < r*r

            p1x, p1y, p2x, p2y, intersects := SegmentCircleIntersection(s1x, s1y, s2x, s2y, cx, cy, r)

            if intersects {
                solutionCntr++
                //Check if points are on circle
                c1x, c1y := p1x-cx, p1y-cy
                deltaLen1 := (c1x*c1x + c1y*c1y) - r*r
                CheckEpsilon(t, deltaLen1, epsilon, "p1 not on circle")

                c2x, c2y := p2x-cx, p2y-cy
                deltaLen2 := (c2x*c2x + c2y*c2y) - r*r
                CheckEpsilon(t, deltaLen2, epsilon, "p2 not on circle")

                // Check if points are on the line through the line segment
                // "cross product" of vector from a segment point to the point
                // and the vector for the segment should be near zero
                vp1x, vp1y := p1x-s1x, p1y-s1y
                crossProd1 := vp1x*segy - vp1y*segx
                CheckEpsilon(t, crossProd1, epsilon, "p1 not on line ")

                vp2x, vp2y := p2x-s1x, p2y-s1y
                crossProd2 := vp2x*segy - vp2y*segx
                CheckEpsilon(t, crossProd2, epsilon, "p2 not on line ")

                // Check if point is between points s1 and s2 on line
                // This means the sign of the dot prod of the segment vector
                // and point to segment end point vectors are opposite for
                // either end.
                wp1x, wp1y := p1x-s2x, p1y-s2y
                dp1v := vp1x*segx + vp1y*segy
                dp1w := wp1x*segx + wp1y*segy
                if (dp1v < 0 && dp1w < 0) || (dp1v > 0 && dp1w > 0) {
                    t.Error("point not contained in segment ", dp1v, dp1w)
                    t.FailNow()
                }

                wp2x, wp2y := p2x-s2x, p2y-s2y
                dp2v := vp2x*segx + vp2y*segy
                dp2w := wp2x*segx + wp2y*segy
                if (dp2v < 0 && dp2w < 0) || (dp2v > 0 && dp2w > 0) {
                    t.Error("point not contained in segment ", dp2v, dp2w)
                    t.FailNow()
                }

                if s1x == s2x && s2y == s1y { //Only one solution
                    // Test that one end of the segment is withing the radius of the circle
                    // and one is not
                    if seg1Inside && seg2Inside {
                        t.Error("Only one solution but both line segment ends inside")
                        t.FailNow()
                    }
                    if !seg1Inside && !seg2Inside {
                        t.Error("Only one solution but both line segment ends outside")
                        t.FailNow()
                    }

                }
            } else { // No intersection, check if both points outside or inside
                if (seg1Inside && !seg2Inside) || (!seg1Inside && seg2Inside) {
                    t.Error("No solution but only one point in radius of circle")
                    t.FailNow()
                }
            }
        }
    }
    t.Log("Tested ", testCntr, " examples and found ", solutionCntr, " solutions.")
}

下面是测试的输出:

=== RUN   TestSegmentCircleIntersection
--- PASS: TestSegmentCircleIntersection (0.00s)
    geom_test.go:105: Tested  40000  examples and found  7343  solutions.

最后,该方法很容易扩展到射线从一点开始,经过另一点并延伸到无穷远的情况,只需测试t > 0或t < 1,而不是两者都测试。

另一个在c#(部分圆类)。 经过测试,工作就像一个魅力。

public class Circle : IEquatable<Circle>
{
    // ******************************************************************
    // The center of a circle
    private Point _center;
    // The radius of a circle
    private double _radius;

   // ******************************************************************
    /// <summary>
    /// Find all intersections (0, 1, 2) of the circle with a line defined by its 2 points.
    /// Using: http://math.stackexchange.com/questions/228841/how-do-i-calculate-the-intersections-of-a-straight-line-and-a-circle
    /// Note: p is the Center.X and q is Center.Y
    /// </summary>
    /// <param name="linePoint1"></param>
    /// <param name="linePoint2"></param>
    /// <returns></returns>
    public List<Point> GetIntersections(Point linePoint1, Point linePoint2)
    {
        List<Point> intersections = new List<Point>();

        double dx = linePoint2.X - linePoint1.X;

        if (dx.AboutEquals(0)) // Straight vertical line
        {
            if (linePoint1.X.AboutEquals(Center.X - Radius) || linePoint1.X.AboutEquals(Center.X + Radius))
            {
                Point pt = new Point(linePoint1.X, Center.Y);
                intersections.Add(pt);
            }
            else if (linePoint1.X > Center.X - Radius && linePoint1.X < Center.X + Radius)
            {
                double x = linePoint1.X - Center.X;

                Point pt = new Point(linePoint1.X, Center.Y + Math.Sqrt(Radius * Radius - (x * x)));
                intersections.Add(pt);

                pt = new Point(linePoint1.X, Center.Y - Math.Sqrt(Radius * Radius - (x * x)));
                intersections.Add(pt);
            }

            return intersections;
        }

        // Line function (y = mx + b)
        double dy = linePoint2.Y - linePoint1.Y;
        double m = dy / dx;
        double b = linePoint1.Y - m * linePoint1.X;

        double A = m * m + 1;
        double B = 2 * (m * b - m * _center.Y - Center.X);
        double C = Center.X * Center.X + Center.Y * Center.Y - Radius * Radius - 2 * b * Center.Y + b * b;

        double discriminant = B * B - 4 * A * C;

        if (discriminant < 0)
        {
            return intersections; // there is no intersections
        }

        if (discriminant.AboutEquals(0)) // Tangeante (touch on 1 point only)
        {
            double x = -B / (2 * A);
            double y = m * x + b;

            intersections.Add(new Point(x, y));
        }
        else // Secant (touch on 2 points)
        {
            double x = (-B + Math.Sqrt(discriminant)) / (2 * A);
            double y = m * x + b;
            intersections.Add(new Point(x, y));

            x = (-B - Math.Sqrt(discriminant)) / (2 * A);
            y = m * x + b;
            intersections.Add(new Point(x, y));
        }

        return intersections;
    }

    // ******************************************************************
    // Get the center
    [XmlElement("Center")]
    public Point Center
    {
        get { return _center; }
        set
        {
            _center = value;
        }
    }

    // ******************************************************************
    // Get the radius
    [XmlElement]
    public double Radius
    {
        get { return _radius; }
        set { _radius = value; }
    }

    //// ******************************************************************
    //[XmlArrayItemAttribute("DoublePoint")]
    //public List<Point> Coordinates
    //{
    //    get { return _coordinates; }
    //}

    // ******************************************************************
    // Construct a circle without any specification
    public Circle()
    {
        _center.X = 0;
        _center.Y = 0;
        _radius = 0;
    }

    // ******************************************************************
    // Construct a circle without any specification
    public Circle(double radius)
    {
        _center.X = 0;
        _center.Y = 0;
        _radius = radius;
    }

    // ******************************************************************
    // Construct a circle with the specified circle
    public Circle(Circle circle)
    {
        _center = circle._center;
        _radius = circle._radius;
    }

    // ******************************************************************
    // Construct a circle with the specified center and radius
    public Circle(Point center, double radius)
    {
        _center = center;
        _radius = radius;
    }

    // ******************************************************************
    // Construct a circle based on one point
    public Circle(Point center)
    {
        _center = center;
        _radius = 0;
    }

    // ******************************************************************
    // Construct a circle based on two points
    public Circle(Point p1, Point p2)
    {
        Circle2Points(p1, p2);
    }

要求:

using System;

namespace Mathematic
{
    public static class DoubleExtension
    {
        // ******************************************************************
        // Base on Hans Passant Answer on:
        // http://stackoverflow.com/questions/2411392/double-epsilon-for-equality-greater-than-less-than-less-than-or-equal-to-gre

        /// <summary>
        /// Compare two double taking in account the double precision potential error.
        /// Take care: truncation errors accumulate on calculation. More you do, more you should increase the epsilon.
        public static bool AboutEquals(this double value1, double value2)
        {
            if (double.IsPositiveInfinity(value1))
                return double.IsPositiveInfinity(value2);

            if (double.IsNegativeInfinity(value1))
                return double.IsNegativeInfinity(value2);

            if (double.IsNaN(value1))
                return double.IsNaN(value2);

            double epsilon = Math.Max(Math.Abs(value1), Math.Abs(value2)) * 1E-15;
            return Math.Abs(value1 - value2) <= epsilon;
        }

        // ******************************************************************
        // Base on Hans Passant Answer on:
        // http://stackoverflow.com/questions/2411392/double-epsilon-for-equality-greater-than-less-than-less-than-or-equal-to-gre

        /// <summary>
        /// Compare two double taking in account the double precision potential error.
        /// Take care: truncation errors accumulate on calculation. More you do, more you should increase the epsilon.
        /// You get really better performance when you can determine the contextual epsilon first.
        /// </summary>
        /// <param name="value1"></param>
        /// <param name="value2"></param>
        /// <param name="precalculatedContextualEpsilon"></param>
        /// <returns></returns>
        public static bool AboutEquals(this double value1, double value2, double precalculatedContextualEpsilon)
        {
            if (double.IsPositiveInfinity(value1))
                return double.IsPositiveInfinity(value2);

            if (double.IsNegativeInfinity(value1))
                return double.IsNegativeInfinity(value2);

            if (double.IsNaN(value1))
                return double.IsNaN(value2);

            return Math.Abs(value1 - value2) <= precalculatedContextualEpsilon;
        }

        // ******************************************************************
        public static double GetContextualEpsilon(this double biggestPossibleContextualValue)
        {
            return biggestPossibleContextualValue * 1E-15;
        }

        // ******************************************************************
        /// <summary>
        /// Mathlab equivalent
        /// </summary>
        /// <param name="dividend"></param>
        /// <param name="divisor"></param>
        /// <returns></returns>
        public static double Mod(this double dividend, double divisor)
        {
            return dividend - System.Math.Floor(dividend / divisor) * divisor;
        }

        // ******************************************************************
    }
}