Line circle intersection
Line circle intersection is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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.
- See Wolfram for 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 {
// hack to get rid of negative zero
// compiler treats 0 and -0 as being same
if p.x == 0 {
p.x = 0
}
if p.y == 0 {
p.y = 0
}
return fmt.Sprintf("(%g, %g)", 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, -5) (3, -2)]
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, -2)]
A circle, center (0, 0) with radius 4, and:
a vertical line containing the points (0, -3) and (0, 6) is/are:
[(0, 4) (0, -4)]
a vertical segment starting at (0, -3) and ending at (0, 6) is/are:
[(0, 4)]
A circle, center (4, 2) with radius 5, and:
a line containing the points (6, 3) and (10, 7) is/are:
[(8, 5) (1, -2)]
a segment starting at (7, 4) and ending at (11, 8) is/are:
[(8, 5)]
Perl 6
Extend solution space to 3D. Reference: this SO question and answers <lang perl6>#!/usr/bin/env perl6
sub LineCircularOBJintersection(@P1, @P2, @Centre, \Radius) {
my @d = @P2 »-« @P1 ; # d my @f = @P1 »-« @Centre ; # c
my \a = [+] @d»²; # d dot d my \b = 2 * ([+] @f »*« @d); # 2 * f dot d my \c = ([+] @f»²) - Radius²; # f dot f - r² my \Δ = b²-(4*a*c); # discriminant
if (Δ < 0) {
return [];
} else {
my (\t1,\t2) = (-b - Δ.sqrt)/(2*a), (-b + Δ.sqrt)/(2*a);
if ( t1 ≥ 0 and t1 ≤ 1 ) or ( t2 ≥ 0 and t2 ≤ 1 ) {
return @P1 »+« ( @P2 »-« @P1 ) »*» t1, @P1 »+« ( @P2 »-« @P1 ) »*» t2
} else {
return []
}
}
}
my \DATA = [
[ <-10 11 > , < 10 -9 > , < 3 -5 > , 3 ], [ <-10 11 > , <-11 12> , < 3 -5 > , 3], [ <3 -2 > , <7 -2> , < 3 -5 > , 3], [ <3 -2 > , <7 -2> , < 0 0> , 4], [ <0 -3 > , <0 6> , < 0 0> , 4], [ <6 3 > , <10 7> , < 4 2> , 5], [ <7 4 > , <11 18> , < 4 2> , 5], [ <3.63 2 −2.26 > , <0.77 2 3.46 > , <1 4 0 > , 4]
];
for DATA {
my @solution = LineCircularOBJintersection $_[0] , $_[1] , $_[2], $_[3]; say "For data set: ", $_; say "Solution(s) is/are: ", @solution.Bool ?? @solution !! "None";
}</lang>
- Output:
For data set: [(-10 11) (10 -9) (3 -5) 3] Solution(s) is/are: [(3 -2) (6 -5)] For data set: [(-10 11) (-11 12) (3 -5) 3] Solution(s) is/are: None For data set: [(3 -2) (7 -2) (3 -5) 3] Solution(s) is/are: [(3 -2) (3 -2)] For data set: [(3 -2) (7 -2) (0 0) 4] Solution(s) is/are: [(-3.4641016151377544 -2) (3.4641016151377544 -2)] For data set: [(0 -3) (0 6) (0 0) 4] Solution(s) is/are: [(0 -4) (0 4)] For data set: [(6 3) (10 7) (4 2) 5] Solution(s) is/are: [(1 -2) (8 5)] For data set: [(7 4) (11 18) (4 2) 5] Solution(s) is/are: [(5.030680985703315 -2.892616550038399) (7.459885052032535 5.60959768211387)] For data set: [(3.63 2 −2.26) (0.77 2 3.46) (1 4 0) 4] Solution(s) is/are: [(3.62828568570857 2 -2.2565713714171394) (0.7717143142914296 2 3.456571371417141)]
Phix
<lang Phix>constant epsilon = 1e-14 -- say atom cx, cy, r, x1, y1, x2, y2
function sq(atom x) return x*x end function
function within(atom x, y)
--
-- checks whether a point is within a segment
-- ie: <-------d1------->
-- <--d2--><---d3---> -- within, d2+d3 ~= d1
-- x1,y1^ ^x,y ^x2,y2
-- vs:
-- <-d2->
-- <-----------d3---------> -- not "", d2+d3 > d1
-- ^x,y - and obviously ditto when x,y is (say) out here^
--
-- (obviously only works when x,y is on the same line as x1,y1 to x2,y2)
--
atom d1 := sqrt(sq(x2-x1) + sq(y2-y1)), -- distance between end-points
d2 := sqrt(sq(x -x1) + sq(y -y1)), -- distance from point to one end
d3 := sqrt(sq(x2-x ) + sq(y2-y )), -- distance from point to other end
delta := (d2 + d3) - d1
return abs(delta) < epsilon -- true if delta is less than a small tolerance
end function
function pf(atom x,y)
return sprintf("(%g,%g)",{x,y})
end function
function intersects(bool bSegment) -- -- Returns the intersection points (if any) of a circle, center (cx,cy) with radius r, -- and line containing the points (x1,y1) and (x2,y2) being either infinite or limited -- to the segment drawn between those points. --
sequence res = {}
atom A = y2 - y1, sqA = sq(A),
B = x1 - x2, sqB = sq(B),
C = x2*y1 - x1*y2, sqC = sq(C),
sqr = r*r-cx*cx-cy*cy,
a := sqA + sqB,
b, c
bool bDivA = false
if abs(B)<epsilon then -- B is zero or close to it
b = 2 * (B*C + A*B*cx - sqA*cy)
c = sqC + 2*A*C*cx - sqA*sqr
bDivA = true -- (and later divide by A instead!)
else
b = 2 * (A*C + A*B*cy - sqB*cx)
c = sqC + 2*B*C*cy - sqB*sqr
end if
atom d := b*b - 4*a*c -- discriminant
if d>=0 then -- (-ve means line & circle do not intersect)
d = sqrt(d)
atom ux,uy, vx,vy
if bDivA then
{uy,vy} = sq_div(sq_sub({+d,-d},b),2*a)
{ux,vx} = sq_div(sq_sub(sq_mul(-B,{uy,vy}),C),A)
else
{ux,vx} = sq_div(sq_sub({+d,-d},b),2*a)
{uy,vy} = sq_div(sq_sub(sq_mul(-A,{ux,vx}),C),B)
end if
if not bSegment or within(ux,uy) then
res = append(res,pf(ux,uy))
end if
if d!=0 and (not bSegment or within(vx,vy)) then
res = append(res,pf(vx,vy))
end if
end if
return res
end function
-- cx cy r x1 y1 x2 y2 bSegment constant tests = {{3,-5,3,{{-10,11, 10,-9,false},
{-10,11,-11,12,true},
{ 3,-2, 7,-2,false}}},
{0, 0,4,{{ 0,-3, 0, 6,false},
{ 0,-3, 0, 6,true}}},
{4, 2,5,{{ 6, 3, 10, 7,false},
{ 7, 4, 11, 8,true}}}}
for t=1 to length(tests) do
{cx, cy, r, sequence lines} = tests[t]
string circle = sprintf("Circle at %s radius %d",{pf(cx,cy),r})
for l=1 to length(lines) do
{x1, y1, x2, y2, bool bSegment} = lines[l]
sequence res = intersects(bSegment)
string ls = iff(bSegment?"segment":" line"),
at = iff(length(res)?"intersect at "&join(res," and ")
:"do not intersect")
printf(1,"%s and %s %s to %s %s.\n",{circle,ls,pf(x1,y1),pf(x2,y2),at})
circle = repeat(' ',length(circle))
end for
end for</lang>
- Output:
Circle at (3,-5) radius 3 and line (-10,11) to (10,-9) intersect at (6,-5) and (3,-2).
and segment (-10,11) to (-11,12) do not intersect.
and line (3,-2) to (7,-2) intersect at (3,-2).
Circle at (0,0) radius 4 and line (0,-3) to (0,6) intersect at (0,4) and (0,-4).
and segment (0,-3) to (0,6) intersect at (0,4).
Circle at (4,2) radius 5 and line (6,3) to (10,7) intersect at (8,5) and (1,-2).
and segment (7,4) to (11,8) intersect at (8,5).
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))