Line-circle intersection tests

// Idea via D3 and MathWorld #include <iostream> #include <cmath> #include <vector> #include <limits> using namespace std; struct Point { double x; double y; }; struct Line { Point p1, p2; }; struct Circle { Point p; double radius; }; int sign(double val) { if(val < 0) return -1; return 1; } void print_points(const vector& v) { for(size_t i = 0; i < v.size(); i++) { cout << "P" << (i+1) << "(" << v[i].x << ", " << v[i].y << ")" << endl; } }; bool intersection_point_lies_on_line(Point curpoint, Point p1, Point p2) { // See StackOverflow double dxc = curpoint.x - p1.x; double dyc = curpoint.y - p1.y; double dxl = p2.x - p1.x; double dyl = p2.y - p1.y; double cross = dxc * dyl - dyc * dxl; double eps = pow(2, -10); if(cross > eps) return false; if(abs(dxl) >= abs(dyl)) { return dxl > 0 ? p1.x <= curpoint.x && curpoint.x <= p2.x : p2.x <= curpoint.x && curpoint.x <= p1.x; } else { return dyl > 0 ? p1.y <= curpoint.y && curpoint.y <= p2.y : p2.y <= curpoint.y && curpoint.y <= p1.y; } }; vector get_intersection_points(double dx, double dy, double dr, double D, double discriminant) { vector v; double dr2 = dr*dr; double discr_sqrt = sqrt(discriminant); double x_num = sign(dy) * dx * discr_sqrt; double y_num = abs(dy) * discr_sqrt; double x1 = (D*dy + x_num) / dr2; double x2 = (D*dy - x_num) / dr2; double y1 = (-D*dx + y_num) / dr2; double y2 = (-D*dx - y_num) / dr2; Point p1 = {x1, y1}; Point p2 = {x2, y2}; v.push_back(p1); v.push_back(p2); return v; }; void line_and_circle_intersect(const Line& line, const Circle& circle) { Point p1 = line.p1; Point p2 = line.p2; double dx = p2.x - p1.x; double dy = p2.y - p1.y; double r = circle.radius; double dr = sqrt(dx*dx + dy*dy); double D = p1.x * p2.y - p2.x * p1.y; double discriminant = r * r * dr * dr - D * D; if(discriminant != 0) { if(discriminant < 0) {cout << "no intersection" << endl;} if(discriminant > 0) { vector points = get_intersection_points(dx, dy, dr, D, discriminant); bool p1_on_line = intersection_point_lies_on_line(points[0], p1, p2); bool p2_on_line = intersection_point_lies_on_line(points[1], p1, p2); if(p1_on_line || p2_on_line) { cout << "intersection" << endl; if(p1_on_line && p2_on_line) {print_points(points); return;} if(p1_on_line) { vector one_point(points.begin(), points.begin() + 1); print_points(one_point); } if(p2_on_line) { vector one_point(points.begin() + 1, points.end()); print_points(one_point); } } else { cout << "no intersection" << endl; } } } else { cout << "line is tangent to the circle" << endl; } }; int main() { Point p1 = {-20, 15}; Point p2 = {20, 24}; Point center = {0, 0}; // implied vector radiuses = {10, 15, 20, 22, 26}; Circle circle; Line line = {p1, p2}; for(auto radius : radiuses) { circle = {center, radius}; cout << "circle radius = " << radius << ": "; line_and_circle_intersect(line, circle); } /* circle radius = 10: no intersection circle radius = 15: line is tangent to the circle circle radius = 20: intersection P1(13.2288, 15) P2(-13.2288, 15) circle radius = 22: intersection P1(16.0935, 15) P2(-16.0935, 15) circle radius = 26: no intersection */ /* If we changed p2.y = 24, the output becomes: circle radius = 10: no intersection circle radius = 15: no intersection circle radius = 20: intersection P1(1.84372, 19.9148) P2(-10.1959, 17.2059) circle radius = 22: intersection P1(6.60308, 20.9857) P2(-14.9553, 16.1351) circle radius = 26: intersection P1(13.1138, 22.4506) */ return 0; }

Unfortunately, this code implies that the center of the circle is at (0, 0). If the center has different coordinates, it may not work (at least this is what I saw in my test).

Here are visualizations of the cases:

no intersectiontangent lineintersectionintersectionno intersectionone point intersectio

When the radius r of the circle is small, we have no line-circle intersection. When it gets too big, we call an additional function to check whether the intersection point lies on the line. (See the case where r = 26 to understand why.)

For simplicity (and to easily determine a tangent point), I used a constant y-coordinate for the two points of the line. But if we changed the y-value of the second point, the code seems to work equally well (at least it seems to match visual observations). In the last case, when the circle radius is r = 26 and the y-coordinate of the second point of the line is set to 24, we see that there is only a single intersection point, which our code seems to recognize correctly.