Line Intersecting Maximal Number of Circles (Circle "Stabbing" Problem)

One of the hardest problems in this year's ACM PACNW regional competition was a computational geometry problem asking you to find a line that intersects the most circles, given a set of disjoint circles. You can read the problem J - Janeway here. This type of problems in computational geometry are known as "stabbing" problems because the line "stabs" the most number of circle. This post describes a O(n2logn)O(n^2 \log n) solution, where nn is the number of circles.

Finding Circle Circle Tangents

First, let's focus on finding tangent lines to two disjoint circles. Trust me, this will come in very handy later! Like most geometry and computational geometry problems, vector algebra is going to be very useful in deriving the formulae here, so remember that all points on the plane are treated as vectors in R2\mathbb{R}^2 and the usual vector operations are applicable to them. Also, you can have a look at this JSXGraph applet to get a geometric sense of what we are doing here, although keep in mind that below we derive the tangents purely algebraically, as opposed to the method used in the applet which is a "ruler and compass" construction.

Assume two circles, CC and CC' are given, with centers at o=(ox,oy)\vec{o} = (o_x, o_y) and o=(ox,oy)\vec{o'} = (o'_x, o'_y) and radii rr and r,r', respectively. Assume the two circles are disjoint and external to each other. Let \ell be a line tangent to both CC and CC' at points p\vec{p} and p.\vec{p'}. Let v=po.\vec{v} = \vec{p} - \vec{o}. Then clearly v=r.|\vec{v}| = r. To find the point on the CC' corresponding to this tangent, we simply notice that

v=po=±rrv. \vec{v'} = \vec{p'} - \vec{o'} = \pm \frac{r'}{r}\vec{v}.

Hence we can identify the tangent line \ell with the vector v\vec{v} since it uniquely identifies a point on CC corresponding to the tangent \ell . Note that the choice of v=+rrv\vec{v'} = +\frac{r'}{r}\vec{v} leads to an external tangent line and the choice of v=rrv\vec{v'} = -\frac{r'}{r}\vec{v} gives an internal one. See figures below.

External and internal tangents to two circles and corresponding vectors.

Since pp\vec{p} - \vec{p'} is on \ell then it is orthogonal v.\vec{v}. Building on this, we have pp=o+vov\vec{p} - \vec{p'} = \vec{o} + \vec{v} - \vec{o'} - \vec{v'} this leads to the following equations:

(o+vov)v=0, (\vec{o} + \vec{v} - \vec{o'} - \vec{v'}) \cdot \vec{v} = \vec{0},

with v=r.| \vec{v} | = r.

Assume v=(x,y).\vec{v} = (x, y). Then the choice of v=rrv\vec{v'} = \frac{r'}{r}\vec{v} leads to the following equation in terms of ox,oy,ox,oy,ro_x, o_y, o'_x, o'_y, r and rr' :

(oxox)x+(oyoy)y+(1rr)(x2+y2)=0, (o_x-o'_x)x + (o_y-o'_y)y + (1 - \frac{r'}{r})(x^2+y^2) = 0,

and v=rrv\vec{v'} = -\frac{r'}{r}\vec{v} gives

(oxox)x+(oyoy)y+(1+rr)(x2+y2)=0. (o_x-o'_x)x + (o_y-o'_y)y + (1 + \frac{r'}{r})(x^2+y^2) = 0.

Substituting x2+y2=r2x^2 + y^2 = r^2 in these equations simplifies the above to

(oxox)x+(oyoy)y+r2±rr=0. (o_x-o'_x)x + (o_y-o'_y)y + r^2 \pm r r'= 0.

Let μ=oxox,\mu = o_x-o'_x, η=oyoy,\eta = o_y - o'_y, and δ=r2±rr,\delta = r^2 \pm r r', with the sign in delta depending on whether we are looking for an internal tangent line or external. Then the equation can be written as

μx+ηy+δ=0. \mu x + \eta y + \delta = 0.

First, assume that η̸=0.\eta \ne 0. Then we can solve for y,y, getting y=μxδη.y=\frac{-\mu x - \delta}{\eta}. Substituting this into x2+y2=r2x^2 + y^2 = r^2 and simplifying gives

(μ2+η2)x2+2δμx+δ2r2η2=0, (\mu^2 + \eta^2)x^2 + 2\delta\mu x + \delta^2 - r^2 \eta^2 = 0,

which has solutions

x=δμ±δ2μ2(μ2+η2)(δ2r2η2)μ2+η2, x = \frac{-\delta \mu \pm \sqrt{\delta^2 \mu^2 - (\mu^2 + \eta^2)(\delta^2 - r^2 \eta^2)}}{\mu^2 + \eta^2},

provided μ2+η2̸=0,\mu^2 + \eta^2 \ne 0, which is granted if the circles are not concentric. Now, once we have xx we can find y=μxδη.y=\frac{-\mu x - \delta}{\eta}.

If η=0\eta = 0 then oy=oyo_y = o'_y and we have symmetry along the xx -axis. This means that

x=δμμ2, x = \frac{-\delta \mu}{\mu^2},

and for each xx we get y=±r2x2.y = \pm \sqrt{r^2 - x^2}.

The above derivation leads directly to the code below to find the tangent lines.

Finding the Line Intersecting Maximal Circles

Now that we know how to find tangents to two circles, let's reduce the circle stabbing problem to the following simpler problem: given a set of disjoint circles C={C1,,Cn}\mathfrak{C} = \{C_1, \ldots, C_n\} and a fixed circle CCC \in \mathfrak{C} find a line \ell that intersects a maximal number of circles in C\mathfrak{C} and is tangent to C.C.

For this reduction, it suffices to show that a maximal line can be found that is tangent to at least one circle. Starting with any maximal line \ell we can find a line \ell that intersects the same number of circles as \ell , but is tangent to at least one circle. Let n\vec{n} be the normal vector to \ell and let ϵ=+ϵn,\ell_\epsilon = \ell + \epsilon \vec{n}, that is \ell shifted ϵ\epsilon units in the direction of its normal. Then for some ϵ\epsilon we must have ϵ\ell_\epsilon tangent to one of the circles it intersects while still intersecting the same number of circles. Informally, this is because as we increase ϵ\epsilon slowly, starting from ϵ=0,\epsilon = 0, eventually ϵ\ell_\epsilon will intersect fewer circles than .\ell. We can pick ϵ\epsilon to be immediately before the first drop in the number of circles intersected.

More formally, let IC()I_\mathfrak{C}(\ell) be the number of circles in C\mathfrak{C} that \ell intersects. Then we simply let

ϵ=sup{ϵ:IC(+ϵn)=IC()}. \epsilon = \sup\{\epsilon' : I_\mathfrak{C}(\ell + \epsilon' \vec{n}) = I_\mathfrak{C}(\ell) \}.

The above ϵ\epsilon is guaranteed to work, for otherwise it is easy to show a contradiction happens.

Given this reduction, we can focus only on lines that are tangent to at least one circle. This gives the following strategy: for each CCC \in \mathfrak{C} find a line C\ell_C that is tangent to CC and intersects a maximal number of circles in C,\mathfrak{C}, and find the globally maximal \ell among all the C.\ell_C. I show that finding the line C\ell_C as described here can be achieved in time O(nlogn)O(n \log n) and hence finding \ell can be achieved in O(n2logn)O(n^2 \log n) time.

To find C,\ell_C, we observe that given another disjoint circle C,C', there will be exactly four line tangent to both of them, two internal, and two external. Suppose that an external tangent line to both CC and CC' goes through point pEp^E on C.C. To this external tangent line corresponds an internal one, on the same side of the line connecting the two centers of the circles. Assume the internal tangent intersects CC at point pI.p^I. Then for any point pp on CC that is between pEp_E and pI,p_I, (here between has to be determined to be clockwise or counter-clockwise based on which side line connect o\vec{o} to o\vec{o'} the two tangents are) the line tangent to CC at pp intersects CC' as well. Hence we can associate with each pair of corresponding external and internal tangent lines an interval [pE,pI].[p^E, p^I].

For any CkCC_k \in \mathfrak{C} and Ck̸=CC_k \ne C calculate intervals [pkE,pkI][p_k^E, p_k^I] and [qkI,qkE],[q_k^I, q_k^E], corresponding to the two external/internal tangent line pairs, paying careful attention to the correct order in which they appear in the interval based on a counter-clockwise inclusion, and add both of them to a collection of intervals Q.Q.

The problem is now reduced to finding a point pp on CC that is in a maximal number of intervals in Q.Q. Note that QQ has exactly 2(n1)2(n-1) intervals in it.

To find the point p,p, first for each interval [p,q]Q,[p, q] \in Q, add the pairs (Angle(p),1)(Angle(p), 1) and (Angle(q),1)(Angle(q), -1) to an array A.A. Here Angle(p)Angle(p) is to denote the angle the point pp makes on the circle C,C, going counter-clockwise and with the right-most point of CC having angle 0,0, and Angle(q)[0,2π).Angle(q) \in [0, 2\pi).

Now, sort the elements in AA based on their first component, and initialize a variable cc with

c1+{(p,q)Q:Angle(q)>Angle(p)}. c \leftarrow 1 + |\{ (p, q) \in Q : Angle(q) > Angle(p) \}|.

We now start at angle θ=0\theta = 0 and simulate moving along CC counter-clockwise and keeping track of how many circles the tangent line at angle θ\theta intersects as we go along. To do this, we use the elements of AA as the events, and at each event aAa \in A we simply set cc+a[1]c \leftarrow c + a[1] where a[1]a[1] is the second component of aa which is equal to 11 or 1-1 depending on whether the corresponding angle is the beginning of where the tangent lines start intersecting another circle, or the end. We then just need to keep track of the angle θ\theta where cc achieves a maximal value, and voilà, we have C\ell_C !

Below is a C++ implementation of the above algorithm. Note that in the code, I ensure that the first pair of external/internal tangents are to the left of the line connecting the centers of CC and CC' and the second pair to its right. This is needed to ensure the start and end angles are computed correctly, and is done by considering the sign of μ\mu in the first case and the sign of η\eta in the second.

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>

using namespace std;

double EPS = 1e-6;

class Point {
public:
    double x, y;
    Point(double x=0.0, double y=0.0) : x(x), y(y) {}
    double angle() const {
        double a = atan2(y, x);
        if (a < 0) {
            a = atan(1) * 8.0 + a;
        }
        return a;
    }
};

class Event {
public:
    double angle;
    double count;
    Event(double angle = 0, int count = 1) : angle(angle), count(count) {}
    bool operator<(const Event &o) const {
        return angle < o.angle;
    }
};

struct CircleCircleTangents {
public:
    Point external[2];
    Point internal[2];
};

class Circle {
public:
    Point center;
    double radius;
    Circle(double x=0.0, double y=0.0, double r=1.0) : radius(r), center(x,y) {}

    // external[0] and internal[0] are guaranteed to be on the left-side of
    // the directed line connecting C1.center to C2.center
    CircleCircleTangents commonTangents(const Circle& C2) const {
        const Circle& C1 = *this;
        double mu = C1.center.x - C2.center.x;
        double eta = C1.center.y - C2.center.y;
        double r1 = C1.radius;
        double r2 = C2.radius;
        double r1r1 = r1 * r1;
        double r1r2 = r1 * r2;
        double delta1 = r1r1 - r1r2;
        double delta2 = r1r1 + r1r2;
        double D = eta*eta + mu*mu;
        CircleCircleTangents result;
        if (abs(eta) < EPS) {
            // Do not divide by eta! Use x^2 + y^2 = r^2 to find y.
            double dmu = mu < 0? -1 : 1;
            double x = (-delta1 * mu) / D;
            double y = -dmu * sqrt(r1r1 - x * x);
            result.external[0].x = x;
            result.external[0].y = y;
            result.external[1].x = x;
            result.external[1].y = -y;
            x = (-delta2 * mu) / D;
            y = -dmu * sqrt(r1r1 - x * x);
            result.internal[0].x = x;
            result.internal[0].y = y;
            result.internal[1].x = x;
            result.internal[1].y = -y;
        } else {
            // Dividing by eta is ok. Use mu*x + eta*y + delta = 0 to find y.
            double mumu = mu * mu;
            double etaeta = eta * eta;
            double dd1 = delta1 * delta1;
            double dd2 = delta2 * delta2;
            double deta = eta < 0? -1 : 1;
            double Delta1 = deta * sqrt(dd1 * mumu - D*(dd1 - etaeta * r1r1));
            double Delta2 = deta * sqrt(dd2 * mumu - D*(dd2 - etaeta * r1r1));
            double x = (-delta1 * mu + Delta1) / D;
            result.external[0].x = x;
            result.external[0].y = -(mu*x + delta1)/eta;
            x = (-delta1 * mu - Delta1) / D;
            result.external[1].x = x;
            result.external[1].y = -(mu*x + delta1)/eta;
            x = (-delta2 * mu + Delta2) / D;
            result.internal[0].x = x;
            result.internal[0].y = -(mu*x + delta2)/eta;
            x = (-delta2 * mu - Delta2) / D;
            result.internal[1].x = x;
            result.internal[1].y = -(mu*x + delta2)/eta;
        }
        return result;
    }
};

bool add_events(vector<Event>& A, const Point& p, const Point& q) {
    double start = p.angle();
    double end = q.angle();
    A.push_back(Event(start, 1));
    A.push_back(Event(end, -1));
    return start > end;
}

// Given a list of circles, returns (m, c, p) where m is the maximal number of
// circles in C any line can intersect, and p is a point on a circle c in C
// such that the tangent line to c at p intersects m circles in C.
int max_intersecting_line(const Circle* C, int n) {
    int global_max = 1;
    vector<Event> A;
    for(int i = 0; i < n; i++) {
        const Circle& c1 = C[i];
        A.clear();
        int local_max = 1;
        for(int j = 0; j < n; j++) {
            if(j == i) continue;
            const Circle& c2 = C[j];
            CircleCircleTangents Q = c1.commonTangents(c2);
            bool t1 = add_events(A, Q.internal[0], Q.external[0]);
            bool t2 = add_events(A, Q.external[1], Q.internal[1]);
            if(t1 || t2) {
                local_max++;
            }
        }

        if (local_max > global_max) {
            global_max = local_max;
        }

        sort(A.begin(), A.end());
        for(int i = 0; i < A.size(); i++) {
            local_max += A[i].count;
            if(local_max > global_max) {
                global_max = local_max;
            }
        }
    }
    return global_max;
}

int main() {
    Circle C[2000];
    int T;
    cin >> T;
    for (int t = 0; t < T; t++) {
        int n;
        cin >> n;
        for (int i = 0; i < n; i++) {
            cin >> C[i].center.x >> C[i].center.y >> C[i].radius;
        }

        cout << max_intersecting_line(C, n) << endl;
    }
    return 0;
}

Comments