Line circle intersection
In plane geometry, a line (or segment) may intersect a circle at 0, 1 or 2 points.
- Task
Implement a method (function, procedure etc.) in your language which takes as parameters:
- the starting point of a line;
- the point where the line ends;
- the center point of a circle;
- the circle's radius; and
- whether the line is a segment or extends to infinity beyond the above points.
The method should return the intersection points (if any) of the circle and the line.
Illustrate your method with some examples (or use the Go examples below).
- References
- See Math Stack Exchange for development of the formulae needed.
Go
<lang go>package main
import (
"fmt" "math"
)
const eps = 1e-14 // say
type point struct{ x, y float64 }
func (p point) String() string {
return fmt.Sprintf("(%f, %f)", p.x, p.y)
}
func sq(x float64) float64 { return x * x }
// Returns the intersection points (if any) of a circle, center 'cp' with radius 'r', // and either an infinite line containing the points 'p1' and 'p2' // or a segment drawn between those points. func intersects(p1, p2, cp point, r float64, segment bool) []point {
var res []point x0, y0 := cp.x, cp.y x1, y1 := p1.x, p1.y x2, y2 := p2.x, p2.y A := y2 - y1 B := x1 - x2 C := x2*y1 - x1*y2 a := sq(A) + sq(B) var b, c float64 var bnz = true if math.Abs(B) >= eps { // if B isn't zero or close to it b = 2 * (A*C + A*B*y0 - sq(B)*x0) c = sq(C) + 2*B*C*y0 - sq(B)*(sq(r)-sq(x0)-sq(y0)) } else { b = 2 * (B*C + A*B*x0 - sq(A)*y0) c = sq(C) + 2*A*C*x0 - sq(A)*(sq(r)-sq(x0)-sq(y0)) bnz = false } d := sq(b) - 4*a*c // discriminant if d < 0 { // line & circle don't intersect return res }
// checks whether a point is within a segment within := func(x, y float64) bool { d1 := math.Sqrt(sq(x2-x1) + sq(y2-y1)) // distance between end-points d2 := math.Sqrt(sq(x-x1) + sq(y-y1)) // distance from point to one end d3 := math.Sqrt(sq(x2-x) + sq(y2-y)) // distance from point to other end delta := d1 - d2 - d3 return math.Abs(delta) < eps // true if delta is less than a small tolerance }
var x, y float64 fx := func() float64 { return -(A*x + C) / B } fy := func() float64 { return -(B*y + C) / A } rxy := func() { if !segment || within(x, y) { res = append(res, point{x, y}) } }
if d == 0 { // line is tangent to circle, so just one intersect at most if bnz { x = -b / (2 * a) y = fx() rxy() } else { y = -b / (2 * a) x = fy() rxy() } } else { // two intersects at most d = math.Sqrt(d) if bnz { x = (-b + d) / (2 * a) y = fx() rxy() x = (-b - d) / (2 * a) y = fx() rxy() } else { y = (-b + d) / (2 * a) x = fy() rxy() y = (-b - d) / (2 * a) x = fy() rxy() } } return res
}
func main() {
cp := point{3, -5} r := 3.0 fmt.Println("The intersection points (if any) between:") fmt.Println("\n A circle, center (3, -5) with radius 3, and:") fmt.Println("\n a line containing the points (-10, 11) and (10, 9) is/are:") fmt.Println(" ", intersects(point{-10, 11}, point{10, -9}, cp, r, false)) fmt.Println("\n a segment starting at (10, -11) and ending at (-11, 12) is/are") fmt.Println(" ", intersects(point{-10, 11}, point{-11, 12}, cp, r, true)) fmt.Println("\n a horizontal line containing the points (3, -2) and (7, -2) is/are:") fmt.Println(" ", intersects(point{3, -2}, point{7, -2}, cp, r, false)) cp = point{0, 0} r = 4.0 fmt.Println("\n A circle, center (0, 0) with radius 4, and:") fmt.Println("\n a vertical line containing the points (0, -3) and (0, 6) is/are:") fmt.Println(" ", intersects(point{0, -3}, point{0, 6}, cp, r, false)) fmt.Println("\n a vertical segment starting at (0, -3) and ending at (0, 6) is/are:") fmt.Println(" ", intersects(point{0, -3}, point{0, 6}, cp, r, true)) cp = point{4, 2} r = 5.0 fmt.Println("\n A circle, center (4, 2) with radius 5, and:") fmt.Println("\n a line containing the points (6, 3) and (10, 7) is/are:") fmt.Println(" ", intersects(point{6, 3}, point{10, 7}, cp, r, false)) fmt.Println("\n a segment starting at (7, 4) and ending at (11, 8) is/are:") fmt.Println(" ", intersects(point{7, 4}, point{11, 8}, cp, r, true))
}</lang>
- Output:
The intersection points (if any) between: A circle, center (3, -5) with radius 3, and: a line containing the points (-10, 11) and (10, 9) is/are: [(6.000000, -5.000000) (3.000000, -2.000000)] a segment starting at (10, -11) and ending at (-11, 12) is/are [] a horizontal line containing the points (3, -2) and (7, -2) is/are: [(3.000000, -2.000000)] A circle, center (0, 0) with radius 4, and: a vertical line containing the points (0, -3) and (0, 6) is/are: [(-0.000000, 4.000000) (0.000000, -4.000000)] a vertical segment starting at (0, -3) and ending at (0, 6) is/are: [(-0.000000, 4.000000)] A circle, center (4, 2) with radius 5, and: a line containing the points (6, 3) and (10, 7) is/are: [(8.000000, 5.000000) (1.000000, -2.000000)] a segment starting at (7, 4) and ending at (11, 8) is/are: [(8.000000, 5.000000)]
zkl
<lang zkl>const EPS=1e-14; // a close-ness to zero
// p1,p2 are (x,y), circle is ( (x,y),r )
fcn intersectLineCircle(p1,p2, circle, segment=False) // assume line {
cx,cy := circle[0].apply("toFloat"); r := circle[1].toFloat(); x1,y1 := p1.apply("toFloat"); x2,y2 := p2.apply("toFloat"); A,B,C,a := (y2 - y1), (x1 - x2), (x2*y1 - x1*y2), (A*A + B*B); b,c,bnz := 0,0,True; if(B.closeTo(0,EPS)){ // B is zero or close to it b = 2.0 * (B*C + A*B*cx - A*A*cy); c = C*C + 2.0*A*C*cx - A*A*(r*r - cx*cx - cy*cy); bnz = False }else{ b = 2.0*( A*C + A*B*cy - B*B*cx ); c = C*C + 2.0*B*C*cy - B*B*( r*r - cx*cx - cy*cy ); } d := b*b - 4.0*a*c; // discriminant if(d<0.0){ // no real solution? zero --> one solution if (d>-0.005) d=0.0; // close enough to zero else return(T); // no intersection } d=d.sqrt();
reg ux,uy, vx,vy; if(bnz){ ux,vx = (-b + d) / (2*a), (-b - d) / (2*a); uy,vy = -(A*ux + C) / B, -(A*vx + C) / B; }else{ uy,vy = (-b + d) / (2*a), (-b - d) / (2*a); ux,vx = -(B*uy + C) / A, -(B*vy + C) / A; }
if(segment){ within:='wrap(x,y){ // is (x,y) within segment p1 p2?
d1:=( (x2 - x1).pow(2) + (y2 - y1).pow(2) ).sqrt(); d2:=( (x - x1).pow(2) + (y - y1).pow(2) ).sqrt(); d3:=( (x2 - x) .pow(2) + (y2 - y) .pow(2) ).sqrt(); (d1 - d2 - d3).closeTo(0,EPS);
};
i1,i2 := within(ux,uy), within(vx,vy); if(d==0) return(if(i1) T(ux,uy) else T); return(T( i1 and T(ux,uy), i2 and T(vx,vy) ).filter()) }
if(d==0) return( T( T(ux,uy) ) ); return( T(ux,uy), T(vx,vy) )
}</lang> <lang zkl>circle:=T( T(3,-5),3 ); p1,p2 := T(-10,11), T( 10,-9); println("Circle @ ",circle); lcpp(p1,p2,circle); p2:=T(-11,12); lcpp(p1,p2,circle,True); p1,p2 := T(3,-2), T(7,-2); lcpp(p1,p2,circle);
circle:=T( T(0,0),4 ); p1,p2 := T(0,-3), T(0,6); println("\nCircle @ ",circle); lcpp(p1,p2,circle); lcpp(p1,p2,circle,True);
circle:=T( T(4,2),5 ); p1,p2 := T(6,3), T(10,7); println("\nCircle @ ",circle); lcpp(p1,p2,circle); p1,p2 := T(7,4), T(11,8); lcpp(p1,p2,circle,True);
fcn lcpp(p1,p2,circle,segment=False){
println(" %s %s -- %s intersects at %s" .fmt(segment and "Segment" or "Line ", p1,p2,intersectLineCircle(p1,p2, circle,segment)));
}</lang>
- Output:
Circle @ L(L(3,-5),3) Line L(-10,11) -- L(10,-9) intersects at L(L(6,-5),L(3,-2)) Segment L(-10,11) -- L(-11,12) intersects at L() Line L(3,-2) -- L(7,-2) intersects at L(L(3,-2)) Circle @ L(L(0,0),4) Line L(0,-3) -- L(0,6) intersects at L(L(0,4),L(0,-4)) Segment L(0,-3) -- L(0,6) intersects at L(L(0,4)) Circle @ L(L(4,2),5) Line L(6,3) -- L(10,7) intersects at L(L(8,5),L(1,-2)) Segment L(7,4) -- L(11,8) intersects at L(L(8,5))