Line circle intersection: Difference between revisions
m (→{{header|C#}}) |
|||
Line 2,149: | Line 2,149: | ||
[(8.0000, 5.0000)] |
[(8.0000, 5.0000)] |
||
</pre> |
</pre> |
||
=={{header|Visual Basic .NET}}== |
|||
{{trans|C#}} |
|||
<lang vbnet>Module Module1 |
|||
Structure Point |
|||
Implements IComparable(Of Point) |
|||
Public Sub New(mx As Double, my As Double) |
|||
X = mx |
|||
Y = my |
|||
End Sub |
|||
Public ReadOnly Property X As Double |
|||
Public ReadOnly Property Y As Double |
|||
Public Function CompareTo(other As Point) As Integer Implements IComparable(Of Point).CompareTo |
|||
Dim c = X.CompareTo(other.X) |
|||
If c <> 0 Then |
|||
Return c |
|||
End If |
|||
Return Y.CompareTo(other.Y) |
|||
End Function |
|||
Public Overrides Function ToString() As String |
|||
Return String.Format("({0}, {1})", X, Y) |
|||
End Function |
|||
End Structure |
|||
Structure Line |
|||
Public Sub New(mp1 As Point, mp2 As Point, Optional segment As Boolean = False) |
|||
If P2.CompareTo(P1) < 0 Then |
|||
P1 = mp2 |
|||
P2 = mp1 |
|||
Else |
|||
P1 = mp1 |
|||
P2 = mp2 |
|||
End If |
|||
IsSegment = segment |
|||
If P1.X = P2.X Then |
|||
Slope = Double.PositiveInfinity |
|||
YIntercept = Double.NaN |
|||
Else |
|||
Slope = (P2.Y - P1.Y) / (P2.X - P1.X) |
|||
YIntercept = P2.Y - Slope * P2.X |
|||
End If |
|||
End Sub |
|||
Public ReadOnly Property P1 As Point |
|||
Public ReadOnly Property P2 As Point |
|||
Public ReadOnly Property Slope As Double |
|||
Public ReadOnly Property YIntercept As Double |
|||
Public ReadOnly Property IsSegment As Boolean |
|||
Public Function IsVertical() As Boolean |
|||
Return P1.X = P2.X |
|||
End Function |
|||
Public Overrides Function ToString() As String |
|||
Return String.Format("[{0}, {1}]", P1, P2) |
|||
End Function |
|||
End Structure |
|||
Structure Circle |
|||
Public Sub New(c As Point, r As Double) |
|||
Center = c |
|||
Radius = r |
|||
End Sub |
|||
Public ReadOnly Property Center As Point |
|||
Public ReadOnly Property Radius As Double |
|||
Public Function X() As Double |
|||
Return Center.X |
|||
End Function |
|||
Public Function Y() As Double |
|||
Return Center.Y |
|||
End Function |
|||
Public Overrides Function ToString() As String |
|||
Return String.Format("{{ C:{0}, R:{1} }}", Center, Radius) |
|||
End Function |
|||
End Structure |
|||
Function Intersection(oc As Circle, ol As Line) As IEnumerable(Of Point) |
|||
Dim LineIntersection = Iterator Function(ic As Circle, il As Line) As IEnumerable(Of Point) |
|||
Dim m = il.Slope |
|||
Dim c = il.YIntercept |
|||
Dim p = ic.X |
|||
Dim q = ic.Y |
|||
Dim r = ic.Radius |
|||
If il.IsVertical Then |
|||
Dim x = il.P1.X |
|||
Dim B = -2 * q |
|||
Dim CC = p * p + q * q - r * r + x * x - 2 * p * x |
|||
Dim D = B * B - 4 * CC |
|||
If D = 0 Then |
|||
Yield New Point(x, -q) |
|||
ElseIf D > 0 Then |
|||
D = Math.Sqrt(D) |
|||
Yield New Point(x, (-B - D) / 2) |
|||
Yield New Point(x, (-B + D) / 2) |
|||
End If |
|||
Else |
|||
Dim A = m * m + 1 |
|||
Dim B = 2 * (m * c - m * q - p) |
|||
Dim CC = p * p + q * q - r * r + c * c - 2 * c * q |
|||
Dim D = B * B - 4 * A * CC |
|||
If D = 0 Then |
|||
Dim x = -B / (2 * A) |
|||
Dim y = m * x + c |
|||
Yield New Point(x, y) |
|||
ElseIf D > 0 Then |
|||
D = Math.Sqrt(D) |
|||
Dim x = (-B - D) / (2 * A) |
|||
Dim y = m * x + c |
|||
Yield New Point(x, y) |
|||
x = (-B + D) / (2 * A) |
|||
y = m * x + c |
|||
Yield New Point(x, y) |
|||
End If |
|||
End If |
|||
End Function |
|||
Dim int = LineIntersection(oc, ol) |
|||
If ol.IsSegment Then |
|||
Return int.Where(Function(p) p.CompareTo(ol.P1) >= 0 AndAlso p.CompareTo(ol.P2) <= 0) |
|||
Else |
|||
Return int |
|||
End If |
|||
End Function |
|||
Sub Print(c As Circle, lines() As Line) |
|||
Console.WriteLine("Circle: {0}", c) |
|||
For Each line In lines |
|||
Console.Write(vbTab) |
|||
If line.IsSegment Then |
|||
Console.Write("Segment: ") |
|||
Else |
|||
Console.Write("Line: ") |
|||
End If |
|||
Console.WriteLine(line) |
|||
Dim points = Intersection(c, line).ToList |
|||
Console.Write(vbTab + vbTab) |
|||
If points.Count = 0 Then |
|||
Console.WriteLine("do not intersect") |
|||
Else |
|||
Console.WriteLine("intersect at {0}", String.Join(" and ", points)) |
|||
End If |
|||
Next |
|||
Console.WriteLine() |
|||
End Sub |
|||
Sub Main() |
|||
Dim c = New Circle(New Point(3, -5), 3) |
|||
Dim lines() As Line = { |
|||
New Line(New Point(-10, 11), New Point(10, -9)), |
|||
New Line(New Point(-10, 11), New Point(-11, 12), True), |
|||
New Line(New Point(3, -2), New Point(7, -2)) |
|||
} |
|||
Print(c, lines) |
|||
c = New Circle(New Point(0, 0), 4) |
|||
lines = { |
|||
New Line(New Point(0, -3), New Point(0, 6)), |
|||
New Line(New Point(0, -3), New Point(0, 6), True) |
|||
} |
|||
Print(c, lines) |
|||
c = New Circle(New Point(4, 2), 5) |
|||
lines = { |
|||
New Line(New Point(6, 3), New Point(10, 7)), |
|||
New Line(New Point(7, 4), New Point(11, 8), True) |
|||
} |
|||
Print(c, lines) |
|||
End Sub |
|||
End Module</lang> |
|||
{{out}} |
|||
<pre>Circle: { C:(3, -5), R:3 } |
|||
Line: [(-10, 11), (10, -9)] |
|||
intersect at (3, -2) and (6, -5) |
|||
Segment: [(-10, 11), (-11, 12)] |
|||
do not intersect |
|||
Line: [(3, -2), (7, -2)] |
|||
intersect at (3, -2) |
|||
Circle: { C:(0, 0), R:4 } |
|||
Line: [(0, -3), (0, 6)] |
|||
intersect at (0, -4) and (0, 4) |
|||
Segment: [(0, -3), (0, 6)] |
|||
intersect at (0, 4) |
|||
Circle: { C:(4, 2), R:5 } |
|||
Line: [(6, 3), (10, 7)] |
|||
intersect at (1, -2) and (8, 5) |
|||
Segment: [(7, 4), (11, 8)] |
|||
intersect at (8, 5)</pre> |
|||
=={{header|zkl}}== |
=={{header|zkl}}== |
Revision as of 01:49, 3 February 2021
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.
C
<lang c>#include <math.h>
- include <stdbool.h>
- include <stdio.h>
const double eps = 1e-14;
typedef struct point_t {
double x, y;
} point;
point make_point(double x, double y) {
point p = { x, y }; return p;
}
void print_point(point p) {
double x = p.x; double y = p.y; if (x == 0) { x = 0; } if (y == 0) { y = 0; } printf("(%g, %g)", x, y);
}
double sq(double x) {
return x * x;
}
bool within(double x1, double y1, double x2, double y2, double x, double y) {
double d1 = sqrt(sq(x2 - x1) + sq(y2 - y1)); // distance between end-points double d2 = sqrt(sq(x - x1) + sq(y - y1)); // distance from point to one end double d3 = sqrt(sq(x2 - x) + sq(y2 - y)); // distance from point to other end double delta = d1 - d2 - d3; return fabs(delta) < eps; // true if delta is less than a small tolerance
}
int rxy(double x1, double y1, double x2, double y2, double x, double y, bool segment) {
if (!segment || within(x1, y1, x2, y2, x, y)) { print_point(make_point(x, y)); return 1; } else { return 0; }
}
double fx(double A, double B, double C, double x) {
return -(A * x + C) / B;
}
double fy(double A, double B, double C, double y) {
return -(B * y + C) / A;
}
// Prints 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. void intersects(point p1, point p2, point cp, double r, bool segment) {
double x0 = cp.x, y0 = cp.y; double x1 = p1.x, y1 = p1.y; double x2 = p2.x, y2 = p2.y; double A = y2 - y1; double B = x1 - x2; double C = x2 * y1 - x1 * y2; double a = sq(A) + sq(B); double b, c, d; bool bnz = true; int cnt = 0;
if (fabs(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 printf("[]\n"); return; }
if (d == 0) { // line is tangent to circle, so just one intersect at most if (bnz) { double x = -b / (2 * a); double y = fx(A, B, C, x); cnt = rxy(x1, y1, x2, y2, x, y, segment); } else { double y = -b / (2 * a); double x = fy(A, B, C, y); cnt = rxy(x1, y1, x2, y2, x, y, segment); } } else { // two intersects at most d = sqrt(d); if (bnz) { double x = (-b + d) / (2 * a); double y = fx(A, B, C, x); cnt = rxy(x1, y1, x2, y2, x, y, segment);
x = (-b - d) / (2 * a); y = fx(A, B, C, x); cnt += rxy(x1, y1, x2, y2, x, y, segment); } else { double y = (-b + d) / (2 * a); double x = fy(A, B, C, y); cnt = rxy(x1, y1, x2, y2, x, y, segment);
y = (-b - d) / (2 * a); x = fy(A, B, C, y); cnt += rxy(x1, y1, x2, y2, x, y, segment); } }
if (cnt <= 0) { printf("[]"); }
}
int main() {
point cp = make_point(3, -5); double r = 3.0; printf("The intersection points (if any) between:\n"); printf(" A circle, center (3, -5) with radius 3, and:\n"); printf(" a line containing the points (-10, 11) and (10, -9) is/are:\n"); printf(" "); intersects(make_point(-10, 11), make_point(10, -9), cp, r, false); printf("\n a segment starting at (-10, 11) and ending at (-11, 12) is/are\n"); printf(" "); intersects(make_point(-10, 11), make_point(-11, 12), cp, r, true); printf("\n a horizontal line containing the points (3, -2) and (7, -2) is/are:\n"); printf(" "); intersects(make_point(3, -2), make_point(7, -2), cp, r, false); printf("\n");
cp = make_point(0, 0); r = 4.0; printf(" A circle, center (0, 0) with radius 4, and:\n"); printf(" a vertical line containing the points (0, -3) and (0, 6) is/are:\n"); printf(" "); intersects(make_point(0, -3), make_point(0, 6), cp, r, false); printf("\n a vertical segment starting at (0, -3) and ending at (0, 6) is/are:\n"); printf(" "); intersects(make_point(0, -3), make_point(0, 6), cp, r, true); printf("\n");
cp = make_point(4,2); r = 5.0; printf(" A circle, center (4, 2) with radius 5, and:\n"); printf(" a line containing the points (6, 3) and (10, 7) is/are:\n"); printf(" "); intersects(make_point(6, 3), make_point(10, 7), cp, r, false); printf("\n a segment starting at (7, 4) and ending at (11, 8) is/are:\n"); printf(" "); intersects(make_point(7, 4), make_point(11, 8), cp, r, true); printf("\n");
return 0;
}</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)
C++
<lang cpp>#include <iostream>
- include <utility>
- include <vector>
using Point = std::pair<double, double>; constexpr auto eps = 1e-14;
std::ostream &operator<<(std::ostream &os, const Point &p) {
auto x = p.first; if (x == 0.0) { x = 0.0; } auto y = p.second; if (y == 0.0) { y = 0.0; } return os << '(' << x << ", " << y << ')';
}
template <typename T> std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
auto itr = v.cbegin(); auto end = v.cend();
os << '['; if (itr != end) { os << *itr; itr = std::next(itr); } while (itr != end) { os << ", " << *itr; itr = std::next(itr); } return os << ']';
}
double sq(double x) {
return x * x;
}
std::vector<Point> intersects(const Point &p1, const Point &p2, const Point &cp, double r, bool segment) {
std::vector<Point> res; auto x0 = cp.first; auto y0 = cp.second; auto x1 = p1.first; auto y1 = p1.second; auto x2 = p2.first; auto y2 = p2.second; auto A = y2 - y1; auto B = x1 - x2; auto C = x2 * y1 - x1 * y2; auto a = sq(A) + sq(B); double b, c; bool bnz = true; if (abs(B) >= eps) { 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; } auto d = sq(b) - 4 * a * c; // discriminant if (d < 0) { return res; }
// checks whether a point is within a segment auto within = [x1, y1, x2, y2](double x, double y) { auto d1 = sqrt(sq(x2 - x1) + sq(y2 - y1)); // distance between end-points auto d2 = sqrt(sq(x - x1) + sq(y - y1)); // distance from point to one end auto d3 = sqrt(sq(x2 - x) + sq(y2 - y)); // distance from point to other end auto delta = d1 - d2 - d3; return abs(delta) < eps; // true if delta is less than a small tolerance };
auto fx = [A, B, C](double x) { return -(A * x + C) / B; };
auto fy = [A, B, C](double y) { return -(B * y + C) / A; };
auto rxy = [segment, &res, within](double x, double y) { if (!segment || within(x, y)) { res.push_back(std::make_pair(x, y)); } };
double x, y; if (d == 0.0) { // line is tangent to circle, so just one intersect at most if (bnz) { x = -b / (2 * a); y = fx(x); rxy(x, y); } else { y = -b / (2 * a); x = fy(y); rxy(x, y); } } else { // two intersects at most d = sqrt(d); if (bnz) { x = (-b + d) / (2 * a); y = fx(x); rxy(x, y); x = (-b - d) / (2 * a); y = fx(x); rxy(x, y); } else { y = (-b + d) / (2 * a); x = fy(y); rxy(x, y); y = (-b - d) / (2 * a); x = fy(y); rxy(x, y); } }
return res;
}
int main() {
std::cout << "The intersection points (if any) between:\n";
auto cp = std::make_pair(3.0, -5.0); auto r = 3.0; std::cout << " A circle, center " << cp << " with radius " << r << ", and:\n";
auto p1 = std::make_pair(-10.0, 11.0); auto p2 = std::make_pair(10.0, -9.0); std::cout << " a line containing the points " << p1 << " and " << p2 << " is/are:\n"; std::cout << " " << intersects(p1, p2, cp, r, false) << '\n';
p2 = std::make_pair(-10.0, 12.0); std::cout << " a segment starting at " << p1 << " and ending at " << p2 << " is/are:\n"; std::cout << " " << intersects(p1, p2, cp, r, true) << '\n';
p1 = std::make_pair(3.0, -2.0); p2 = std::make_pair(7.0, -2.0); std::cout << " a horizontal line containing the points " << p1 << " and " << p2 << " is/are:\n"; std::cout << " " << intersects(p1, p2, cp, r, false) << '\n';
cp = std::make_pair(0.0, 0.0); r = 4.0; std::cout << " A circle, center " << cp << " with radius " << r << ", and:\n";
p1 = std::make_pair(0.0, -3.0); p2 = std::make_pair(0.0, 6.0); std::cout << " a vertical line containing the points " << p1 << " and " << p2 << " is/are:\n"; std::cout << " " << intersects(p1, p2, cp, r, false) << '\n'; std::cout << " a vertical segment containing the points " << p1 << " and " << p2 << " is/are:\n"; std::cout << " " << intersects(p1, p2, cp, r, true) << '\n';
cp = std::make_pair(4.0, 2.0); r = 5.0; std::cout << " A circle, center " << cp << " with radius " << r << ", and:\n";
p1 = std::make_pair(6.0, 3.0); p2 = std::make_pair(10.0, 7.0); std::cout << " a line containing the points " << p1 << " and " << p2 << " is/are:\n"; std::cout << " " << intersects(p1, p2, cp, r, false) << '\n';
p1 = std::make_pair(7.0, 4.0); p2 = std::make_pair(11.0, 8.0); std::cout << " a segment starting at " << p1 << " and ending at " << p2 << " is/are:\n"; std::cout << " " << intersects(p1, p2, cp, r, true) << '\n';
return 0;
}</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 (-10, 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 containing the points (0, -3) and (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)]
C#
<lang Csharp>using System; using System.Collections.Generic; using System.Linq;
public class Program {
public static void Main() { Circle circle = ((3, -5), 3); Line[] lines = { ((-10, 11), (10, -9)), ((-10, 11), (-11, 12), true), ((3, -2), (7, -2)) }; Print(circle, lines); circle = ((0, 0), 4); lines = new Line[] { ((0, -3), (0, 6)), ((0, -3), (0, 6), true) }; Print(circle, lines); circle = ((4, 2), 5); lines = new Line[] { ((6, 3), (10, 7)), ((7, 4), (11, 8), true) }; Print(circle, lines); } static void Print(Circle circle, Line[] lines) { Console.WriteLine($"Circle: {circle}"); foreach (var line in lines) { Console.WriteLine($"\t{(line.IsSegment ? "Segment:" : "Line:")} {line}"); var points = Intersection(circle, line).ToList(); Console.WriteLine(points.Count == 0 ? "\t\tdo not intersect" : "\t\tintersect at " + string.Join(" and ", points)); } Console.WriteLine(); } static IEnumerable<Point> Intersection(Circle circle, Line line) { var intersection = LineIntersection(circle, line); return line.IsSegment ? intersection.Where(p => p.CompareTo(line.P1) >= 0 && p.CompareTo(line.P2) <= 0) : intersection;
static IEnumerable<Point> LineIntersection(Circle circle, Line line) { double x, y, A, B, C, D; var (m, c) = (line.Slope, line.YIntercept); var (p, q, r) = (circle.X, circle.Y, circle.Radius);
if (line.IsVertical) { x = line.P1.X; B = -2 * q; C = p * p + q * q - r * r + x * x - 2 * p * x; D = B * B - 4 * C; if (D == 0) yield return (x, -q); else if (D > 0) { D = Math.Sqrt(D); yield return (x, (-B - D) / 2); yield return (x, (-B + D) / 2); } } else { A = m * m + 1; B = 2 * (m * c - m * q - p); C = p * p + q * q - r * r + c * c - 2 * c * q; D = B * B - 4 * A * C; if (D == 0) { x = -B / (2 * A); y = m * x + c; yield return (x, y); } else if (D > 0) { D = Math.Sqrt(D); x = (-B - D) / (2 * A); y = m * x + c; yield return (x, y); x = (-B + D) / (2 * A); y = m * x + c; yield return (x, y); } } }
} readonly struct Point : IComparable<Point> { public Point(double x, double y) => (X, Y) = (x, y); public static implicit operator Point((double x, double y) p) => new Point(p.x, p.y); public double X { get; } public double Y { get; } public int CompareTo(Point other) { int c = X.CompareTo(other.X); if (c != 0) return c; return Y.CompareTo(other.Y); } public override string ToString() => $"({X}, {Y})"; } readonly struct Line { public Line(Point p1, Point p2, bool isSegment = false) { (P1, P2) = p2.CompareTo(p1) < 0 ? (p2, p1) : (p1, p2); IsSegment = isSegment; if (p1.X == p2.X) (Slope, YIntercept) = (double.PositiveInfinity, double.NaN); else { Slope = (P2.Y - P1.Y) / (P2.X - P1.X); YIntercept = P2.Y - Slope * P2.X; } } public static implicit operator Line((Point p1, Point p2) l) => new Line(l.p1, l.p2); public static implicit operator Line((Point p1, Point p2, bool isSegment) l) => new Line(l.p1, l.p2, l.isSegment); public Point P1 { get; } public Point P2 { get; } public double Slope { get; } public double YIntercept { get; } public bool IsSegment { get; } public bool IsVertical => P1.X == P2.X; public override string ToString() => $"[{P1}, {P2}]"; } readonly struct Circle { public Circle(Point center, double radius) => (Center, Radius) = (center, radius); public static implicit operator Circle((Point center, double radius) c) => new Circle(c.center, c.radius); public Point Center { get; } public double Radius { get; } public double X => Center.X; public double Y => Center.Y; public override string ToString() => $"{{ C:{Center}, R:{Radius} }}"; }
}</lang>
- Output:
Circle: { C:(3, -5), R:3 } Line: [(-10, 11), (10, -9)] intersect at (3, -2) and (6, -5) Segment: [(-11, 12), (-10, 11)] do not intersect Line: [(3, -2), (7, -2)] intersect at (3, -2) Circle: { C:(0, 0), R:4 } Line: [(0, -3), (0, 6)] intersect at (0, -4) and (0, 4) Segment: [(0, -3), (0, 6)] intersect at (0, 4) Circle: { C:(4, 2), R:5 } Line: [(6, 3), (10, 7)] intersect at (1, -2) and (8, 5) Segment: [(7, 4), (11, 8)] intersect at (8, 5)
D
<lang d>import std.format; import std.math; import std.stdio;
immutable EPS = 1e-14;
struct Point {
private double x; private double y;
public this(double x, double y) { this.x = x; this.y = y; }
public double getX() { return x; }
public double getY() { return y; }
void toString(scope void delegate(const(char)[]) sink, FormatSpec!char fmt) const { double mx = x; double my = y;
// eliminate negative zero if (mx == 0.0) { mx = 0.0; }
// eliminate negative zero if (my == 0.0) { my = 0.0; }
sink("("); formatValue(sink, mx, fmt); sink(", "); formatValue(sink, my, fmt); sink(")"); }
}
auto sq(T)(T x) {
return x * x;
}
auto intersects(const Point p1, const Point p2, const Point cp, double r, bool segment) {
auto x0 = cp.x; auto y0 = cp.y; auto x1 = p1.x; auto y1 = p1.y; auto x2 = p2.x; auto y2 = p2.y;
auto A = y2 - y1; auto B = x1 - x2; auto C = x2 * y1 - x1 * y2;
auto a = sq(A) + sq(B); double b, c;
bool bnz = true;
Point[] res;
if (abs(B) >= EPS) { 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; }
auto d = sq(b) - 4 * a * c; // discriminant if (d < 0) { return res; }
auto within(double x, double y) { auto d1 = sqrt(sq(x2 - x1) + sq(y2 - y1)); // distance between end-points auto d2 = sqrt(sq(x - x1) + sq(y - y1)); // distance from point to one end auto d3 = sqrt(sq(x2 - x) + sq(y2 - y)); // distance from point to other end auto delta = d1 - d2 - d3; return abs(delta) < EPS; // true if delta is less than a small tolerance }
auto fx(double x) { return -(A * x + C) / B; }
auto fy(double y) { return -(B * y + C) / A; }
auto rxy(double x, double y) { if (!segment || within(x, y)) { res ~= Point(x, y); } }
double x, y; if (d == 0.0) { // line is tangent to circle, so just one intersect at most if (bnz) { x = -b / (2 * a); y = fx(x); rxy(x, y); } else { y = -b / (2 * a); x = fy(y); rxy(x, y); } } else { // two intersects at most d = sqrt(d); if (bnz) { x = (-b + d) / (2 * a); y = fx(x); rxy(x, y); x = (-b - d) / (2 * a); y = fx(x); rxy(x, y); } else { y = (-b + d) / (2 * a); x = fy(y); rxy(x, y); y = (-b - d) / (2 * a); x = fy(y); rxy(x, y); } }
return res;
}
void main() {
writeln("The intersection points (if any) between:");
auto cp = Point(3.0, -5.0); auto r = 3.0; writeln(" A circle, center ", cp, " with radius ", r, ", and:");
auto p1 = Point(-10.0, 11.0); auto p2 = Point(10.0, -9.0); writeln(" a line containing the points ", p1, " and ", p2, " is/are:"); writeln(" ", intersects(p1, p2, cp, r, false));
p2 = Point(-10.0, 12.0); writeln(" a segment starting at ", p1, " and ending at ", p2, " is/are:"); writeln(" ", intersects(p1, p2, cp, r, true));
p1 = Point(3.0, -2.0); p2 = Point(7.0, -2.0); writeln(" a horizontal line containing the points ", p1, " and ", p2, " is/are:"); writeln(" ", intersects(p1, p2, cp, r, false));
cp = Point(0.0, 0.0); r = 4.0; writeln(" A circle, center ", cp, " with radius ", r, ", and:");
p1 = Point(0.0, -3.0); p2 = Point(0.0, 6.0); writeln(" a vertical line containing the points ", p1, " and ", p2, " is/are:"); writeln(" ", intersects(p1, p2, cp, r, false)); writeln(" a vertical segment containing the points ", p1, " and ", p2, " is/are:"); writeln(" ", intersects(p1, p2, cp, r, true));
cp = Point(4.0, 2.0); r = 5.0; writeln(" A circle, center ", cp, " with radius ", r, ", and:");
p1 = Point(6.0, 3.0); p2 = Point(10.0, 7.0); writeln(" a line containing the points ", p1, " and ", p2, " is/are:"); writeln(" ", intersects(p1, p2, cp, r, false));
p1 = Point(7.0, 4.0); p2 = Point(11.0, 8.0); writeln(" a segment starting at ", p1, " and ending at ", p2, " is/are:"); writeln(" ", intersects(p1, p2, 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 (-10, 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 containing the points (0, -3) and (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)]
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)]
Haskell
<lang haskell>import Data.Tuple.Curry
main :: IO () main =
mapM_ putStrLn $ concatMap (("" :) . uncurryN task) [ ((-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)) ]
task :: (Double, Double)
-> (Double, Double) -> ((Double, Double), Double) -> [String]
task pt1 pt2 circle@(pt3@(a3, b3), r) = [line, segment]
where xs = map fun $ lineCircleIntersection pt1 pt2 circle ys = map fun $ segmentCircleIntersection pt1 pt2 circle to x = (fromIntegral . round $ 100 * x) / 100 fun (x, y) = (to x, to y) yo = show . fun start = "Intersection: Circle " ++ yo pt3 ++ " " ++ show (to r) ++ " and " end = yo pt1 ++ " " ++ yo pt2 ++ ": " line = start ++ "Line " ++ end ++ show xs segment = start ++ "Segment " ++ end ++ show ys
segmentCircleIntersection
:: (Double, Double) -> (Double, Double) -> ((Double, Double), Double) -> [(Double, Double)]
segmentCircleIntersection pt1 pt2 circle =
filter (go p1 p2) $ lineCircleIntersection pt1 pt2 circle where [p1, p2] | pt1 < pt2 = [pt1, pt2] | otherwise = [pt2, pt1] go (x, y) (u, v) (i, j) | x == u = y <= j && j <= v | otherwise = x <= i && i <= u
lineCircleIntersection
:: (Double, Double) -> (Double, Double) -> ((Double, Double), Double) -> [(Double, Double)]
lineCircleIntersection (a1, b1) (a2, b2) ((a3, b3), r) = go delta
where (x1, x2) = (a1 - a3, a2 - a3) (y1, y2) = (b1 - b3, b2 - b3) (dx, dy) = (x2 - x1, y2 - y1) drdr = dx * dx + dy * dy d = x1 * y2 - x2 * y1 delta = r * r * drdr - d * d sqrtDelta = sqrt delta (sgnDy, absDy) = (sgn dy, abs dy) u1 = (d * dy + sgnDy * dx * sqrtDelta) / drdr u2 = (d * dy - sgnDy * dx * sqrtDelta) / drdr v1 = (-d * dx + absDy * sqrtDelta) / drdr v2 = (-d * dx - absDy * sqrtDelta) / drdr go x | 0 > x = [] | 0 == x = [(u1 + a3, v1 + b3)] | otherwise = [(u1 + a3, v1 + b3), (u2 + a3, v2 + b3)]
sgn :: Double -> Double sgn x
| 0 > x = -1 | otherwise = 1</lang>
- Output:
Intersection: Circle (3.0,-5.0) 3.0 and Line (-10.0,11.0) (10.0,-9.0): [(3.0,-2.0),(6.0,-5.0)] Intersection: Circle (3.0,-5.0) 3.0 and Segment (-10.0,11.0) (10.0,-9.0): [(3.0,-2.0),(6.0,-5.0)] Intersection: Circle (3.0,-5.0) 3.0 and Line (-10.0,11.0) (-11.0,12.0): [(3.0,-2.0),(6.0,-5.0)] Intersection: Circle (3.0,-5.0) 3.0 and Segment (-10.0,11.0) (-11.0,12.0): [] Intersection: Circle (3.0,-5.0) 3.0 and Line (3.0,-2.0) (7.0,-2.0): [(3.0,-2.0)] Intersection: Circle (3.0,-5.0) 3.0 and Segment (3.0,-2.0) (7.0,-2.0): [(3.0,-2.0)] Intersection: Circle (0.0,0.0) 4.0 and Line (3.0,-2.0) (7.0,-2.0): [(3.46,-2.0),(-3.46,-2.0)] Intersection: Circle (0.0,0.0) 4.0 and Segment (3.0,-2.0) (7.0,-2.0): [(3.46,-2.0)] Intersection: Circle (0.0,0.0) 4.0 and Line (0.0,-3.0) (0.0,6.0): [(0.0,4.0),(0.0,-4.0)] Intersection: Circle (0.0,0.0) 4.0 and Segment (0.0,-3.0) (0.0,6.0): [(0.0,4.0)] Intersection: Circle (4.0,2.0) 5.0 and Line (6.0,3.0) (10.0,7.0): [(8.0,5.0),(1.0,-2.0)] Intersection: Circle (4.0,2.0) 5.0 and Segment (6.0,3.0) (10.0,7.0): [(8.0,5.0)] Intersection: Circle (4.0,2.0) 5.0 and Line (7.0,4.0) (11.0,18.0): [(7.46,5.61),(5.03,-2.89)] Intersection: Circle (4.0,2.0) 5.0 and Segment (7.0,4.0) (11.0,18.0): [(7.46,5.61)]
Java
<lang java>import java.util.*; import java.awt.geom.*;
public class LineCircleIntersection {
public static void main(String[] args) { try { demo(); } catch (Exception e) { e.printStackTrace(); } }
private static void demo() throws NoninvertibleTransformException { Point2D center = makePoint(3, -5); double radius = 3.0; System.out.println("The intersection points (if any) between:"); System.out.println("\n A circle, center (3, -5) with radius 3, and:"); System.out.println("\n a line containing the points (-10, 11) and (10, -9) is/are:"); System.out.println(" " + toString(intersection(makePoint(-10, 11), makePoint(10, -9), center, radius, false))); System.out.println("\n a segment starting at (-10, 11) and ending at (-11, 12) is/are"); System.out.println(" " + toString(intersection(makePoint(-10, 11), makePoint(-11, 12), center, radius, true))); System.out.println("\n a horizontal line containing the points (3, -2) and (7, -2) is/are:"); System.out.println(" " + toString(intersection(makePoint(3, -2), makePoint(7, -2), center, radius, false))); center.setLocation(0, 0); radius = 4.0; System.out.println("\n A circle, center (0, 0) with radius 4, and:"); System.out.println("\n a vertical line containing the points (0, -3) and (0, 6) is/are:"); System.out.println(" " + toString(intersection(makePoint(0, -3), makePoint(0, 6), center, radius, false))); System.out.println("\n a vertical segment starting at (0, -3) and ending at (0, 6) is/are:"); System.out.println(" " + toString(intersection(makePoint(0, -3), makePoint(0, 6), center, radius, true))); center.setLocation(4, 2); radius = 5.0; System.out.println("\n A circle, center (4, 2) with radius 5, and:"); System.out.println("\n a line containing the points (6, 3) and (10, 7) is/are:"); System.out.println(" " + toString(intersection(makePoint(6, 3), makePoint(10, 7), center, radius, false))); System.out.println("\n a segment starting at (7, 4) and ending at (11, 8) is/are:"); System.out.println(" " + toString(intersection(makePoint(7, 4), makePoint(11, 8), center, radius, true))); }
private static Point2D makePoint(double x, double y) { return new Point2D.Double(x, y); }
// // If center of the circle is at the origin and the line is horizontal, // it's easy to calculate the points of intersection, so to handle the // general case, we convert the input to a coordinate system where the // center of the circle is at the origin and the line is horizontal, // then convert the points of intersection back to the original // coordinate system. // public static List<Point2D> intersection(Point2D p1, Point2D p2, Point2D center, double radius, boolean isSegment) throws NoninvertibleTransformException { List<Point2D> result = new ArrayList<>(); double dx = p2.getX() - p1.getX(); double dy = p2.getY() - p1.getY(); AffineTransform trans = AffineTransform.getRotateInstance(dx, dy); trans.invert(); trans.translate(-center.getX(), -center.getY()); Point2D p1a = trans.transform(p1, null); Point2D p2a = trans.transform(p2, null); double y = p1a.getY(); double minX = Math.min(p1a.getX(), p2a.getX()); double maxX = Math.max(p1a.getX(), p2a.getX()); if (y == radius || y == -radius) { if (!isSegment || (0 <= maxX && 0 >= minX)) { p1a.setLocation(0, y); trans.inverseTransform(p1a, p1a); result.add(p1a); } } else if (y < radius && y > -radius) { double x = Math.sqrt(radius * radius - y * y); if (!isSegment || (-x <= maxX && -x >= minX)) { p1a.setLocation(-x, y); trans.inverseTransform(p1a, p1a); result.add(p1a); } if (!isSegment || (x <= maxX && x >= minX)) { p2a.setLocation(x, y); trans.inverseTransform(p2a, p2a); result.add(p2a); } } return result; }
public static String toString(Point2D point) { return String.format("(%g, %g)", point.getX(), point.getY()); }
public static String toString(List<Point2D> points) { StringBuilder str = new StringBuilder("["); for (int i = 0, n = points.size(); i < n; ++i) { if (i > 0) str.append(", "); str.append(toString(points.get(i))); } str.append("]"); return str.toString(); }
}</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: [(3.00000, -2.00000), (6.00000, -5.00000)] 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.00000, -2.00000)] A circle, center (0, 0) with radius 4, and: a vertical line containing the points (0, -3) and (0, 6) is/are: [(0.00000, -4.00000), (0.00000, 4.00000)] a vertical segment starting at (0, -3) and ending at (0, 6) is/are: [(0.00000, 4.00000)] A circle, center (4, 2) with radius 5, and: a line containing the points (6, 3) and (10, 7) is/are: [(1.00000, -2.00000), (8.00000, 5.00000)] a segment starting at (7, 4) and ending at (11, 8) is/are: [(8.00000, 5.00000)]
Julia
Uses the circles and points from the Go example. <lang julia>using Luxor
const centers = [Point(3, -5), Point(0, 0), Point(4, 2)] const rads = [3, 4, 5] const lins = [
[Point(-10, 11), Point(10, -9)], [Point(-10, 11), Point(-11, 12)], [Point(3, -2), Point(7, -2)], [Point(0, -3), Point(0, 6)], [Point(6, 3), Point(10, 7)], [Point(7, 4), Point(11, 8)],
]
println("Center", " "^9, "Radius", " "^4, "Line P1", " "^14, "Line P2", " "^7,
"Segment? Intersect 1 Intersect 2")
for (cr, l, extended) in [(1, 1, true), (1, 2, false), (1, 3, false),
(2, 4, true), (2, 4, false), (3, 5, true), (3, 6, false)] tup = intersectionlinecircle(lins[l][1], lins[l][2], centers[cr], rads[cr]) v = [p for p in tup[2:end] if extended || ispointonline(p, lins[l][1], lins[l][2])] println(rpad(centers[cr], 17), rads[cr], " "^3, rpad(lins[l][1], 21), rpad(lins[l][2], 19), rpad(!extended, 8), isempty(v) ? "" : length(v) == 2 ? rpad(v[1], 18) * string(v[2]) : v[1])
end
</lang>
- Output:
Center Radius Line P1 Line P2 Segment? Intersect 1 Intersect 2 Point(3.0, -5.0) 3 Point(-10.0, 11.0) Point(10.0, -9.0) false Point(6.0, -5.0) Point(3.0, -2.0) Point(3.0, -5.0) 3 Point(-10.0, 11.0) Point(-11.0, 12.0) true Point(3.0, -5.0) 3 Point(3.0, -2.0) Point(7.0, -2.0) true Point(3.0, -2.0) Point(0.0, 0.0) 4 Point(0.0, -3.0) Point(0.0, 6.0) false Point(0.0, 4.0) Point(0.0, -4.0) Point(0.0, 0.0) 4 Point(0.0, -3.0) Point(0.0, 6.0) true Point(0.0, 4.0) Point(4.0, 2.0) 5 Point(6.0, 3.0) Point(10.0, 7.0) false Point(8.0, 5.0) Point(1.0, -2.0) Point(4.0, 2.0) 5 Point(7.0, 4.0) Point(11.0, 8.0) true Point(8.0, 5.0)
Kotlin
<lang scala>import kotlin.math.absoluteValue import kotlin.math.sqrt
const val eps = 1e-14
class Point(val x: Double, val y: Double) {
override fun toString(): String { var xv = x if (xv == 0.0) { xv = 0.0 } var yv = y if (yv == 0.0) { yv = 0.0 } return "($xv, $yv)" }
}
fun sq(x: Double): Double {
return x * x
}
fun intersects(p1: Point, p2: Point, cp: Point, r: Double, segment: Boolean): MutableList<Point> {
val res = mutableListOf<Point>() val x0 = cp.x val y0 = cp.y val x1 = p1.x val y1 = p1.y val x2 = p2.x val y2 = p2.y val A = y2 - y1 val B = x1 - x2 val C = x2 * y1 - x1 * y2 val a = sq(A) + sq(B) val b: Double val c: Double var bnz = true if (B.absoluteValue >= eps) { 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 } var d = sq(b) - 4 * a * c // discriminant if (d < 0) { return res }
// checks whether a point is within a segment fun within(x: Double, y: Double): Boolean { val d1 = sqrt(sq(x2 - x1) + sq(y2 - y1)) // distance between end-points val d2 = sqrt(sq(x - x1) + sq(y - y1)) // distance from point to one end val d3 = sqrt(sq(x2 - x) + sq(y2 - y)) // distance from point to other end val delta = d1 - d2 - d3 return delta.absoluteValue < eps // true if delta is less than a small tolerance }
var x = 0.0 fun fx(): Double { return -(A * x + C) / B }
var y = 0.0 fun fy(): Double { return -(B * y + C) / A }
fun rxy() { if (!segment || within(x, y)) { res.add(Point(x, y)) } }
if (d == 0.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 = 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
}
fun main() {
println("The intersection points (if any) between:")
var cp = Point(3.0, -5.0) var r = 3.0 println(" A circle, center $cp with radius $r, and:")
var p1 = Point(-10.0, 11.0) var p2 = Point(10.0, -9.0) println(" a line containing the points $p1 and $p2 is/are:") println(" ${intersects(p1, p2, cp, r, false)}")
p2 = Point(-10.0, 12.0) println(" a segment starting at $p1 and ending at $p2 is/are:") println(" ${intersects(p1, p2, cp, r, true)}")
p1 = Point(3.0, -2.0) p2 = Point(7.0, -2.0) println(" a horizontal line containing the points $p1 and $p2 is/are:") println(" ${intersects(p1, p2, cp, r, false)}")
cp = Point(0.0, 0.0) r = 4.0 println(" A circle, center $cp with radius $r, and:")
p1 = Point(0.0, -3.0) p2 = Point(0.0, 6.0) println(" a vertical line containing the points $p1 and $p2 is/are:") println(" ${intersects(p1, p2, cp, r, false)}") println(" a vertical segment containing the points $p1 and $p2 is/are:") println(" ${intersects(p1, p2, cp, r, true)}")
cp = Point(4.0, 2.0) r = 5.0 println(" A circle, center $cp with radius $r, and:")
p1 = Point(6.0, 3.0) p2 = Point(10.0, 7.0) println(" a line containing the points $p1 and $p2 is/are:") println(" ${intersects(p1, p2, cp, r, false)}")
p1 = Point(7.0, 4.0) p2 = Point(11.0, 8.0) println(" a segment starting at $p1 and ending at $p2 is/are:") println(" ${intersects(p1, p2, cp, r, true)}")
}</lang>
- Output:
The intersection points (if any) between: A circle, center (3.0, -5.0) with radius 3.0, and: a line containing the points (-10.0, 11.0) and (10.0, -9.0) is/are: [(6.0, -5.0), (3.0, -2.0)] a segment starting at (-10.0, 11.0) and ending at (-10.0, 12.0) is/are: [] a horizontal line containing the points (3.0, -2.0) and (7.0, -2.0) is/are: [(3.0, -2.0)] A circle, center (0.0, 0.0) with radius 4.0, and: a vertical line containing the points (0.0, -3.0) and (0.0, 6.0) is/are: [(0.0, 4.0), (0.0, -4.0)] a vertical segment containing the points (0.0, -3.0) and (0.0, 6.0) is/are: [(0.0, 4.0)] A circle, center (4.0, 2.0) with radius 5.0, and: a line containing the points (6.0, 3.0) and (10.0, 7.0) is/are: [(8.0, 5.0), (1.0, -2.0)] a segment starting at (7.0, 4.0) and ending at (11.0, 8.0) is/are: [(8.0, 5.0)]
Lua
<lang lua>EPS = 1e-14
function pts(p)
local x, y = p.x, p.y if x == 0 then x = 0 end if y == 0 then y = 0 end return "(" .. x .. ", " .. y .. ")"
end
function lts(pl)
local str = "[" for i,p in pairs(pl) do if i > 1 then str = str .. ", " end str = str .. pts(p) end return str .. "]"
end
function sq(x)
return x * x
end
function intersects(p1, p2, cp, r, segment)
local res = {} local x0, y0 = cp.x, cp.y local x1, y1 = p1.x, p1.y local x2, y2 = p2.x, p2.y local A = y2 - y1 local B = x1 - x2 local C = x2 * y1 - x1 * y2 local a = sq(A) + sq(B) local b, c local bnz = true if math.abs(B) >= EPS then 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 end local d = sq(b) - 4 * a * c -- discriminant if d < 0 then return res end
-- checks whether a point is within a segment function within(x, y) local d1 = math.sqrt(sq(x2 - x1) + sq(y2 - y1)) -- distance between end-points local d2 = math.sqrt(sq(x - x1) + sq(y - y1)) -- distance from point to one end local d3 = math.sqrt(sq(x2 - x) + sq(y2 - y)) -- distance from point to other end local delta = d1 - d2 - d3 return math.abs(delta) < EPS end
function fx(x) return -(A * x + C) / B end
function fy(y) return -(B * y + C) / A end
function rxy(x, y) if not segment or within(x, y) then table.insert(res, {x=x, y=y}) end end
local x, y if d == 0 then -- line is tangent to circle, so just one intersect at most if bnz then x = -b / (2 * a) y = fx(x) rxy(x, y) else y = -b / (2 * a) x = fy(y) rxy(x, y) end else -- two intersects at most d = math.sqrt(d) if bnz then x = (-b + d) / (2 * a) y = fx(x) rxy(x, y) x = (-b - d) / (2 * a) y = fx(x) rxy(x, y) else y = (-b + d) / (2 * a) x = fy(y) rxy(x, y) y = (-b - d) / (2 * a) x = fy(y) rxy(x, y) end end
return res
end
function main()
print("The intersection points (if any) between:")
local cp = {x=3, y=-5} local r = 3 print(" A circle, center " .. pts(cp) .. " with radius " .. r .. ", and:")
local p1 = {x=-10, y=11} local p2 = {x=10, y=-9} print(" a line containing the points " .. pts(p1) .. " and " .. pts(p2) .. " is/are:") print(" " .. lts(intersects(p1, p2, cp, r, false)))
p2 = {x=-10, y=12} print(" a segment starting at " .. pts(p1) .. " and ending at " .. pts(p2) .. " is/are:") print(" " .. lts(intersects(p1, p2, cp, r, true)))
p1 = {x=3, y=-2} p2 = {x=7, y=-2} print(" a horizontal line containing the points " .. pts(p1) .. " and " .. pts(p2) .. " is/are:") print(" " .. lts(intersects(p1, p2, cp, r, false)))
cp = {x=0, y=0} r = 4 print(" A circle, center " .. pts(cp) .. " with radius " .. r .. ", and:")
p1 = {x=0, y=-3} p2 = {x=0, y=6} print(" a vertical line containing the points " .. pts(p1) .. " and " .. pts(p2) .. " is/are:") print(" " .. lts(intersects(p1, p2, cp, r, false))) print(" a vertical segment starting at " .. pts(p1) .. " and ending at " .. pts(p2) .. " is/are:") print(" " .. lts(intersects(p1, p2, cp, r, true)))
cp = {x=4, y=2} r = 5 print(" A circle, center " .. pts(cp) .. " with radius " .. r .. ", and:")
p1 = {x=6, y=3} p2 = {x=10, y=7} print(" a line containing the points " .. pts(p1) .. " and " .. pts(p2) .. " is/are:") print(" " .. lts(intersects(p1, p2, cp, r, false)))
p1 = {x=7, y=4} p2 = {x=11, y=8} print(" a segment starting at " .. pts(p1) .. " and ending at " .. pts(p2) .. " is/are:") print(" " .. lts(intersects(p1, p2, cp, r, true)))
end
main()</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 (-10, 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
<lang perl>use strict; use warnings; use feature 'say'; use List::Util 'sum';
sub find_intersection {
my($P1, $P2, $center, $radius) = @_; my @d = ($$P2[0] - $$P1[0], $$P2[1] - $$P1[1]); my @f = ($$P1[0] - $$center[0], $$P1[1] - $$center[1]); my $a = sum map { $_**2 } @d; my $b = 2 * ($f[0]*$d[0] + $f[1]*$d[1]); my $c = sum(map { $_**2 } @f) - $radius**2; my $D = $b**2 - 4*$a*$c;
return unless $D >= 0; my ($t1, $t2) = ( (-$b - sqrt $D) / (2*$a), (-$b + sqrt $D) / (2*$a) ); return unless $t1 >= 0 and $t1 <= 1 or $t2 >= 0 and $t2 <= 1;
my ($dx, $dy) = ($$P2[0] - $$P1[0], $$P2[1] - $$P1[1]); return ([$dx*$t1 + $$P1[0], $dy*$t1 + $$P1[1]], [$dx*$t2 + $$P1[0], $dy*$t2 + $$P1[1]])
}
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 ],
);
sub rnd { map { sprintf('%.2f', $_) =~ s/\.00//r } @_ }
for my $d (@data) {
my @solution = find_intersection @$d[0] , @$d[1] , @$d[2], @$d[3]; say 'For input: ' . join ', ', (map { '('. join(',', @$_) .')' } @$d[0,1,2]), @$d[3]; say 'Solutions: ' . (@solution > 1 ? join ', ', map { '('. join(',', rnd @$_) .')' } @solution : 'None'); say ;
}</lang>
- Output:
For input: (-10,11), (10,-9), (3,-5), 3 Solutions: (3,-2), (6,-5) For input: (-10,11), (-11,12), (3,-5), 3 Solutions: None For input: (3,-2), (7,-2), (3,-5), 3 Solutions: (3,-2), (3,-2) For input: (3,-2), (7,-2), (0,0), 4 Solutions: (-3.46,-2), (3.46,-2) For input: (0,-3), (0,6), (0,0), 4 Solutions: (0,-4), (0,4) For input: (6,3), (10,7), (4,2), 5 Solutions: (1,-2), (8,5) For input: (7,4), (11,18), (4,2), 5 Solutions: (5.03,-2.89), (7.46,5.61)
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).
Raku
(formerly Perl 6) Extend solution space to 3D. Reference: this SO question and answers <lang 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 0 ≤ t1|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 ], [ <5 2 −2.26 >, <0.77 2 4>, <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: [(5 2 −2.26) (0.77 2 4) (1 4 0) 4] Solution(s) is/are: [(4.2615520237084015 2 -1.1671668246843006) (1.13386504516801 2 3.461514141193441)]
REXX
The formulae used for this REXX version were taken from the MathWorld webpage: circle line intersection. <lang rexx>/*REXX program calculates where (or if) a line intersects (or tengents) a cirle. */ /*───────────────────────────────────── line= x1,y1 x2,y2; circle is at 0,0, radius=r*/ parse arg x1 y1 x2 y2 cx cy r . /*obtain optional arguments from the CL*/ if x1== | x1=="," then x1= 0 /*Not specified? Then use the default.*/ if y1== | y1=="," then y1= -3 /* " " " " " " */ if x2== | x2=="," then x2= 0 /* " " " " " " */ if y2== | y2=="," then y2= 6 /* " " " " " " */ if cx== | cx=="," then cx= 0 /* " " " " " " */ if cy== | cy=="," then cy= 0 /* " " " " " " */ if r == | r =="," then r = 4 /* " " " " " " */ x_1= x1; x1= x1 + cx; y_1= y1; y1= y1 + cy x_2= x2; x2= x2 + cx; y_2= y2; y2= y2 + cy dx= x2 - x1; dy= y2 - y1
dr2= dx**2 + dy**2 D= x1 * y2 - x2 * y1; r2= r**2; D2= D**2 $= sqrt(r2 * dr2 - D2)
ix1= ( D * dy + sgn(dy) * dx * $) / dr2 ix2= ( D * dy - sgn(dy) * dx * $) / dr2 iy1= (-D * dx + abs(dy) * $) / dr2 iy2= (-D * dx - abs(dy) * $) / dr2 incidence= (r2 * dr2 - D2) / 1 say 'incidence=' incidence
@potla= 'points on the line are: '
if incidence<0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1 say "The line doesn't intersect the circle with radius: " r end
if incidence=0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1 say "The line is tangent to circle with radius: " r end
if incidence>0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1 say "The line is secant to circle with radius: " r end
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ sgn: procedure; if arg(1)<0 then return -1; return 1 /*──────────────────────────────────────────────────────────────────────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2 do j=0 while h>9; m.j= h; h= h%2 + 1; end /*j*/ do k=j+5 to 0 by -1; numeric digits m.k; g= (g+x/g) *.5; end /*k*/; return g</lang>
- output when using the default inputs:
incidence= 1296 points on the line are: (0,-3) and (0,6) are: 0,4 The line is secant to circle with radius: 4
Ruby
<lang ruby>EPS = 1e-14
def sq(x)
return x * x
end
def intersects(p1, p2, cp, r, segment)
res = [] (x0, y0) = cp (x1, y1) = p1 (x2, y2) = p2 aa = y2 - y1 bb = x1 - x2 cc = x2 * y1 - x1 * y2 a = sq(aa) + sq(bb) if bb.abs >= EPS then b = 2 * (aa * cc + aa * bb * y0 - sq(bb) * x0) c = sq(cc) + 2 * bb * cc * y0 - sq(bb) * (sq(r) - sq(x0) - sq(y0)) bnz = true else b = 2 * (bb * cc + aa * bb * x0 - sq(aa) * y0) c = sq(cc) + 2 * aa * cc * x0 - sq(aa) * (sq(r) - sq(x0) - sq(y0)) bnz = false end d = sq(b) - 4 * a * c # disciminant if d < 0 then return res end
# checks whether a point is within a segment within = ->(x, y) { 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 delta.abs < EPS # true if delta is less than a small tolerance }
fx = ->(x) { return -(aa * x + cc) / bb }
fy = ->(y) { return -(bb * y + cc) / aa }
rxy = ->(x, y) { if not segment or within.call(x, y) then if x == 0.0 then x = 0.0 end if y == 0.0 then y = 0.0 end res << [x, y] end }
if d == 0.0 then # line is tangent to circle, so just one intersect at most if bnz then x = -b / (2 * a) y = fx.call(x) rxy.call(x, y) else y = -b / (2 * a) x = fy.call(y) rxy.call(x, y) end else # two intersects at most d = Math.sqrt(d) if bnz then x = (-b + d) / (2 * a) y = fx.call(x) rxy.call(x, y) x = (-b - d) / (2 * a) y = fx.call(x) rxy.call(x, y) else y = (-b + d) / (2 * a) x = fy.call(y) rxy.call(x, y) y = (-b - d) / (2 * a) x = fy.call(y) rxy.call(x, y) end end
return res
end
def main
print "The intersection points (if any) between:\n"
cp = [3.0, -5.0] r = 3.0 print " A circle, center %s with radius %f, and:\n" % [cp, r]
p1 = [-10.0, 11.0] p2 = [10.0, -9.0] print " a line containing the points %s and %s is/are:\n" % [p1, p2] print " %s\n" % [intersects(p1, p2, cp, r, false)]
p2 = [-10.0, 12.0] print " a segment starting at %s and ending at %s is/are:\n" % [p1, p2] print " %s\n" % [intersects(p1, p2, cp, r, true)]
p1 = [3.0, -2.0] p2 = [7.0, -2.0] print " a horizontal line containing the points %s and %s is/are:\n" % [p1, p2] print " %s\n" % [intersects(p1, p2, cp, r, false)]
cp = [0.0, 0.0] r = 4.0 print " A circle, center %s with radius %f, and:\n" % [cp, r]
p1 = [0.0, -3.0] p2 = [0.0, 6.0] print " a vertical line containing the points %s and %s is/are:\n" % [p1, p2] print " %s\n" % [intersects(p1, p2, cp, r, false)] print " a vertical line segment containing the points %s and %s is/are:\n" % [p1, p2] print " %s\n" % [intersects(p1, p2, cp, r, true)]
cp = [4.0, 2.0] r = 5.0 print " A circle, center %s with radius %f, and:\n" % [cp, r]
p1 = [6.0, 3.0] p2 = [10.0, 7.0] print " a line containing the points %s and %s is/are:\n" % [p1, p2] print " %s\n" % [intersects(p1, p2, cp, r, false)]
p1 = [7.0, 4.0] p2 = [11.0, 8.0] print " a segment starting at %s and ending at %s is/are:\n" % [p1, p2] print " %s\n" % [intersects(p1, p2, cp, r, true)]
end
main()</lang>
- Output:
The intersection points (if any) between: A circle, center [3.0, -5.0] with radius 3.000000, and: a line containing the points [-10.0, 11.0] and [10.0, -9.0] is/are: [[6.0, -5.0], [3.0, -2.0]] a segment starting at [-10.0, 11.0] and ending at [-10.0, 12.0] is/are: [] a horizontal line containing the points [3.0, -2.0] and [7.0, -2.0] is/are: [[3.0, -2.0]] A circle, center [0.0, 0.0] with radius 4.000000, and: a vertical line containing the points [0.0, -3.0] and [0.0, 6.0] is/are: [[0.0, 4.0], [0.0, -4.0]] a vertical line segment containing the points [0.0, -3.0] and [0.0, 6.0] is/are: [[0.0, 4.0]] A circle, center [4.0, 2.0] with radius 5.000000, and: a line containing the points [6.0, 3.0] and [10.0, 7.0] is/are: [[8.0, 5.0], [1.0, -2.0]] a segment starting at [7.0, 4.0] and ending at [11.0, 8.0] is/are: [[8.0, 5.0]]
Swift
<lang swift>import Foundation import CoreGraphics
func lineCircleIntersection(start: NSPoint, end: NSPoint, center: NSPoint,
radius: CGFloat, segment: Bool) -> [NSPoint] { var result: [NSPoint] = [] let angle = atan2(end.y - start.y, end.x - start.x) var at = AffineTransform(rotationByRadians: angle) at.invert() at.translate(x: -center.x, y: -center.y) let p1 = at.transform(start) let p2 = at.transform(end) let minX = min(p1.x, p2.x), maxX = max(p1.x, p2.x) let y = p1.y at.invert() func addPoint(x: CGFloat, y: CGFloat) { if !segment || (x <= maxX && x >= minX) { result.append(at.transform(NSMakePoint(x, y))) } } if y == radius || y == -radius { addPoint(x: 0, y: y) } else if y < radius && y > -radius { let x = (radius * radius - y * y).squareRoot() addPoint(x: -x, y: y) addPoint(x: x, y: y) } return result
}
func toString(points: [NSPoint]) -> String {
var result = "[" result += points.map{String(format: "(%.4f, %.4f)", $0.x, $0.y)}.joined(separator: ", ") result += "]" return result
}
var center = NSMakePoint(3, -5) var radius: CGFloat = 3
print("The intersection points (if any) between:") print("\n A circle, center (3, -5) with radius 3, and:") print("\n a line containing the points (-10, 11) and (10, -9) is/are:") var points = lineCircleIntersection(start: NSMakePoint(-10, 11), end: NSMakePoint(10, -9),
center: center, radius: radius, segment: false)
print(" \(toString(points: points))") print("\n a segment starting at (-10, 11) and ending at (-11, 12) is/are") points = lineCircleIntersection(start: NSMakePoint(-10, 11), end: NSMakePoint(-11, 12),
center: center, radius: radius, segment: true)
print(" \(toString(points: points))") print("\n a horizontal line containing the points (3, -2) and (7, -2) is/are:") points = lineCircleIntersection(start: NSMakePoint(3, -2), end: NSMakePoint(7, -2),
center: center, radius: radius, segment: false)
print(" \(toString(points: points))")
center.x = 0 center.y = 0 radius = 4
print("\n A circle, center (0, 0) with radius 4, and:") print("\n a vertical line containing the points (0, -3) and (0, 6) is/are:") points = lineCircleIntersection(start: NSMakePoint(0, -3), end: NSMakePoint(0, 6),
center: center, radius: radius, segment: false)
print(" \(toString(points: points))") print("\n a vertical segment starting at (0, -3) and ending at (0, 6) is/are:") points = lineCircleIntersection(start: NSMakePoint(0, -3), end: NSMakePoint(0, 6),
center: center, radius: radius, segment: true)
print(" \(toString(points: points))")
center.x = 4 center.y = 2 radius = 5
print("\n A circle, center (4, 2) with radius 5, and:") print("\n a line containing the points (6, 3) and (10, 7) is/are:") points = lineCircleIntersection(start: NSMakePoint(6, 3), end: NSMakePoint(10, 7),
center: center, radius: radius, segment: false)
print(" \(toString(points: points))") print("\n a segment starting at (7, 4) and ending at (11, 8) is/are:") points = lineCircleIntersection(start: NSMakePoint(7, 4), end: NSMakePoint(11, 8),
center: center, radius: radius, segment: true)
print(" \(toString(points: points))")</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: [(3.0000, -2.0000), (6.0000, -5.0000)] 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.0000, -2.0000)] A circle, center (0, 0) with radius 4, and: a vertical line containing the points (0, -3) and (0, 6) is/are: [(-0.0000, -4.0000), (0.0000, 4.0000)] a vertical segment starting at (0, -3) and ending at (0, 6) is/are: [(0.0000, 4.0000)] A circle, center (4, 2) with radius 5, and: a line containing the points (6, 3) and (10, 7) is/are: [(1.0000, -2.0000), (8.0000, 5.0000)] a segment starting at (7, 4) and ending at (11, 8) is/are: [(8.0000, 5.0000)]
Visual Basic .NET
<lang vbnet>Module Module1
Structure Point Implements IComparable(Of Point)
Public Sub New(mx As Double, my As Double) X = mx Y = my End Sub
Public ReadOnly Property X As Double Public ReadOnly Property Y As Double
Public Function CompareTo(other As Point) As Integer Implements IComparable(Of Point).CompareTo Dim c = X.CompareTo(other.X) If c <> 0 Then Return c End If Return Y.CompareTo(other.Y) End Function
Public Overrides Function ToString() As String Return String.Format("({0}, {1})", X, Y) End Function End Structure
Structure Line Public Sub New(mp1 As Point, mp2 As Point, Optional segment As Boolean = False) If P2.CompareTo(P1) < 0 Then P1 = mp2 P2 = mp1 Else P1 = mp1 P2 = mp2 End If IsSegment = segment If P1.X = P2.X Then Slope = Double.PositiveInfinity YIntercept = Double.NaN Else Slope = (P2.Y - P1.Y) / (P2.X - P1.X) YIntercept = P2.Y - Slope * P2.X End If End Sub
Public ReadOnly Property P1 As Point Public ReadOnly Property P2 As Point Public ReadOnly Property Slope As Double Public ReadOnly Property YIntercept As Double Public ReadOnly Property IsSegment As Boolean
Public Function IsVertical() As Boolean Return P1.X = P2.X End Function
Public Overrides Function ToString() As String Return String.Format("[{0}, {1}]", P1, P2) End Function End Structure
Structure Circle Public Sub New(c As Point, r As Double) Center = c Radius = r End Sub
Public ReadOnly Property Center As Point Public ReadOnly Property Radius As Double
Public Function X() As Double Return Center.X End Function
Public Function Y() As Double Return Center.Y End Function
Public Overrides Function ToString() As String Return String.Format("{{ C:{0}, R:{1} }}", Center, Radius) End Function End Structure
Function Intersection(oc As Circle, ol As Line) As IEnumerable(Of Point) Dim LineIntersection = Iterator Function(ic As Circle, il As Line) As IEnumerable(Of Point) Dim m = il.Slope Dim c = il.YIntercept Dim p = ic.X Dim q = ic.Y Dim r = ic.Radius
If il.IsVertical Then Dim x = il.P1.X Dim B = -2 * q Dim CC = p * p + q * q - r * r + x * x - 2 * p * x Dim D = B * B - 4 * CC If D = 0 Then Yield New Point(x, -q) ElseIf D > 0 Then D = Math.Sqrt(D) Yield New Point(x, (-B - D) / 2) Yield New Point(x, (-B + D) / 2) End If Else Dim A = m * m + 1 Dim B = 2 * (m * c - m * q - p) Dim CC = p * p + q * q - r * r + c * c - 2 * c * q Dim D = B * B - 4 * A * CC If D = 0 Then Dim x = -B / (2 * A) Dim y = m * x + c Yield New Point(x, y) ElseIf D > 0 Then D = Math.Sqrt(D) Dim x = (-B - D) / (2 * A) Dim y = m * x + c Yield New Point(x, y) x = (-B + D) / (2 * A) y = m * x + c Yield New Point(x, y) End If End If End Function
Dim int = LineIntersection(oc, ol) If ol.IsSegment Then Return int.Where(Function(p) p.CompareTo(ol.P1) >= 0 AndAlso p.CompareTo(ol.P2) <= 0) Else Return int End If End Function
Sub Print(c As Circle, lines() As Line) Console.WriteLine("Circle: {0}", c) For Each line In lines Console.Write(vbTab) If line.IsSegment Then Console.Write("Segment: ") Else Console.Write("Line: ") End If Console.WriteLine(line)
Dim points = Intersection(c, line).ToList
Console.Write(vbTab + vbTab) If points.Count = 0 Then Console.WriteLine("do not intersect") Else Console.WriteLine("intersect at {0}", String.Join(" and ", points)) End If Next Console.WriteLine() End Sub
Sub Main() Dim c = New Circle(New Point(3, -5), 3) Dim lines() As Line = { New Line(New Point(-10, 11), New Point(10, -9)), New Line(New Point(-10, 11), New Point(-11, 12), True), New Line(New Point(3, -2), New Point(7, -2)) } Print(c, lines)
c = New Circle(New Point(0, 0), 4) lines = { New Line(New Point(0, -3), New Point(0, 6)), New Line(New Point(0, -3), New Point(0, 6), True) } Print(c, lines)
c = New Circle(New Point(4, 2), 5) lines = { New Line(New Point(6, 3), New Point(10, 7)), New Line(New Point(7, 4), New Point(11, 8), True) } Print(c, lines) End Sub
End Module</lang>
- Output:
Circle: { C:(3, -5), R:3 } Line: [(-10, 11), (10, -9)] intersect at (3, -2) and (6, -5) Segment: [(-10, 11), (-11, 12)] do not intersect Line: [(3, -2), (7, -2)] intersect at (3, -2) Circle: { C:(0, 0), R:4 } Line: [(0, -3), (0, 6)] intersect at (0, -4) and (0, 4) Segment: [(0, -3), (0, 6)] intersect at (0, 4) Circle: { C:(4, 2), R:5 } Line: [(6, 3), (10, 7)] intersect at (1, -2) and (8, 5) Segment: [(7, 4), (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))