Пошук спільних дотичних до двох кіл

Дано два кола. Потрібно знайти всі їхні спільні дотичні, тобто всі такі прямі, які торкаються обох кіл одночасно.

Описаний алгоритм працюватиме й у випадку, коли одне (чи обидва) кола вироджуються в точки. Тому цей алгоритм можна застосувати і для пошуку дотичних до кола, що проходять через задану точку.

Коли підходить цей алгоритм?

• Вам потрібні спільні дотичні до двох кіл (або дотичні до одного кола через зовнішню точку — як вироджений випадок із колом нульового радіуса)?
• Ви готові обробити окремо вироджений випадок збіжних кіл (нескінченно багато дотичних), оскільки алгоритм його не покриває?
• Вам потрібні рівняння самих дотичних прямих, а не точки перетину кола з прямою? (якщо потрібен перетин кола й прямої → Перетин кола та прямої)

Кількість спільних дотичних

Кількість спільних дотичних до двох кіл може дорівнювати 0, 1, 2, 3, 4 і бути нескінченною. Подивіться на зображення для різних випадків.

Тут ми не розглядатимемо вироджені випадки, тобто коли кола збігаються (у цьому випадку вони мають нескінченно багато спільних дотичних), або одне коло лежить усередині іншого (у цьому випадку вони не мають спільних дотичних, а якщо кола дотичні, то є одна спільна дотична).

У більшості випадків два кола мають чотири спільні дотичні.

Якщо кола дотикаються, то вони матимуть три спільні дотичні, але це можна розуміти як вироджений випадок: ніби дві дотичні збіглися.

Більше того, описаний нижче алгоритм працюватиме у випадку, коли одне чи обидва кола мають нульовий радіус: у цьому випадку буде, відповідно, дві або одна спільна дотична.

Підсумовуючи, ми завжди шукатимемо чотири дотичні для всіх випадків, окрім випадку з нескінченною кількістю дотичних (випадок нескінченної кількості дотичних потрібно обробляти окремо, і тут він не обговорюється). У вироджених випадках деякі з дотичних збігатимуться, але попри це ці випадки теж вписуються в загальну картину.

Алгоритм

Заради простоти алгоритму будемо вважати, без втрати загальності, що центр першого кола має координати $(0, 0)$. (Якщо це не так, то цього можна досягти простим зсувом усієї картинки, а після знаходження розв'язку — зворотним зсувом отриманих прямих.)

Позначимо через $r_1$ і $r_2$ радіуси першого та другого кіл, а через $(v_x,v_y)$ — координати центра другого кола, точки $v$, відмінної від початку координат. (Зауваження: ми не розглядаємо випадок, коли обидва кола однакові.)

Щоб розв'язати задачу, підійдемо до неї суто алгебраїчно. Нам потрібно знайти всі прямі вигляду $ax + by + c = 0$, які лежать на відстані $r_1$ від початку координат і на відстані $r_2$ від точки $v$. Крім того, накладемо умову нормування прямої: сума квадратів коефіцієнтів має дорівнювати одиниці (це необхідно, інакше одній і тій самій прямій відповідатиме нескінченно багато представлень вигляду $ax + by + c = 0$). Загалом отримуємо таку систему рівнянь для шуканих $a, b, c$:

$$\begin{align} a^2 + b^2 &= 1 \\ \mid a \cdot 0 + b \cdot 0 + c \mid &= r_1 \\ \mid a \cdot v_x + b \cdot v_y + c \mid &= r_2 \end{align}$$

Щоб позбутися модуля, зауважимо, що в цій системі є лише чотири способи розкрити модуль. Усі ці способи можна розглянути в загальному випадку, якщо розуміти розкриття модуля як те, що коефіцієнт у правій частині може бути помножений на -1. Інакше кажучи, переходимо до такої системи:

$$\begin{align} a^2 + b^2 &= 1 \\ c &= \pm r_1 \\ a \cdot v_x + b \cdot v_y + c &= \pm r_2 \end{align}$$

Увівши позначення $d_1 = \pm r_1$ та $d_2 = \pm r_2$, доходимо висновку, що система повинна мати чотири розв'язки:

$$\begin{align} a^2 + b^2 &= 1 \\ c &= d_1 \\ a \cdot v_x + b \cdot v_y + c &= d_2 \end{align}$$

Розв'язування цієї системи зводиться до розв'язування квадратного рівняння. Ми опустимо всі громіздкі обчислення й одразу наведемо готову відповідь:

$$\begin{align} a &= {( d_2 - d_1 ) v_x \pm v_y \sqrt{v_x^2 + v_y^2-(d_2-d_1)^2} \over {v_x^2 + v_y^2} } \\ b &= {( d_2 - d_1 ) v_y \pm v_x \sqrt{v_x^2 + v_y^2-(d_2-d_1)^2} \over {v_x^2 + v_y^2} } \\ c &= d_1 \end{align}$$

Загалом ми отримали вісім розв'язків замість чотирьох. Однак легко зрозуміти, звідки виникають зайві розв'язки: насправді в останній системі достатньо взяти лише один розв'язок (наприклад, перший). Справді, геометричний зміст того, що ми беремо $\pm r_1$ та $\pm r_2$, зрозумілий: ми насправді перебираємо, з якого боку від кожного кола проходить пряма. Тому два способи, що виникають при розв'язуванні останньої системи, є надлишковими: достатньо вибрати один із двох розв'язків (тільки, звісно, в усіх чотирьох випадках потрібно вибирати те саме сімейство розв'язків).

Останнє, що ми ще не розглянули, — це як зсунути прямі у випадку, коли перше коло початково не було розташоване в початку координат. Однак тут усе просто: з лінійності рівняння прямої випливає, що значення $a \cdot x_0 + b \cdot y_0$ (де $x_0$ та $y_0$ — координати початкового центра першого кола) потрібно відняти від коефіцієнта $c$.

Реалізація

Спочатку опишемо всі необхідні структури даних та інші допоміжні означення:

C++

struct pt {
    double x, y;

    pt operator- (pt p) {
        pt res = { x-p.x, y-p.y };
        return res;
    }
};

struct circle : pt {
    double r;
};

struct line {
    double a, b, c;
};

const double EPS = 1E-9;

double sqr (double a) {
    return a * a;
}

Python

from dataclasses import dataclass

EPS = 1e-9


@dataclass
class Pt:
    x: float
    y: float

    def __sub__(self, p: "Pt") -> "Pt":
        return Pt(self.x - p.x, self.y - p.y)


@dataclass
class Circle:
    x: float
    y: float
    r: float


@dataclass
class Line:
    a: float
    b: float
    c: float


def sqr(a: float) -> float:
    return a * a

TypeScript

const EPS = 1e-9;

interface Pt {
    x: number;
    y: number;
}

interface Circle {
    x: number;
    y: number;
    r: number;
}

interface Line {
    a: number;
    b: number;
    c: number;
}

function sub(a: Pt, b: Pt): Pt {
    return { x: a.x - b.x, y: a.y - b.y };
}

function sqr(a: number): number {
    return a * a;
}

Go

package main

import "math"

const EPS = 1e-9

type Pt struct {
    X, Y float64
}

func (p Pt) Sub(q Pt) Pt {
    return Pt{p.X - q.X, p.Y - q.Y}
}

type Circle struct {
    X, Y, R float64
}

type Line struct {
    A, B, C float64
}

func sqr(a float64) float64 {
    return a * a
}

Тоді сам розв'язок можна записати так (де основна функція для виклику — друга, а перша функція є допоміжною):

C++

void tangents (pt c, double r1, double r2, vector<line> & ans) {
    double r = r2 - r1;
    double z = sqr(c.x) + sqr(c.y);
    double d = z - sqr(r);
    if (d < -EPS)  return;
    d = sqrt (abs (d));
    line l;
    l.a = (c.x * r + c.y * d) / z;
    l.b = (c.y * r - c.x * d) / z;
    l.c = r1;
    ans.push_back (l);
}

vector<line> tangents (circle a, circle b) {
    vector<line> ans;
    for (int i=-1; i<=1; i+=2)
        for (int j=-1; j<=1; j+=2)
            tangents (b-a, a.r*i, b.r*j, ans);
    for (size_t i=0; i<ans.size(); ++i)
        ans[i].c -= ans[i].a * a.x + ans[i].b * a.y;
    return ans;
}

Python

import math


def tangents(c: Pt, r1: float, r2: float, ans: list[Line]) -> None:
    r = r2 - r1
    z = sqr(c.x) + sqr(c.y)
    d = z - sqr(r)
    if d < -EPS:
        return
    d = math.sqrt(abs(d))
    ans.append(Line(
        a=(c.x * r + c.y * d) / z,
        b=(c.y * r - c.x * d) / z,
        c=r1,
    ))


def find_tangents(a: Circle, b: Circle) -> list[Line]:
    ans: list[Line] = []
    # Перебираємо знаки ±r1 та ±r2 (з якого боку проходить пряма)
    for i in (-1, 1):
        for j in (-1, 1):
            tangents(Pt(b.x - a.x, b.y - a.y), a.r * i, b.r * j, ans)
    # Зворотний зсув: повертаємо центр першого кола з початку координат
    for l in ans:
        l.c -= l.a * a.x + l.b * a.y
    return ans

TypeScript

function tangents(c: Pt, r1: number, r2: number, ans: Line[]): void {
    const r = r2 - r1;
    const z = sqr(c.x) + sqr(c.y);
    let d = z - sqr(r);
    if (d < -EPS) return;
    d = Math.sqrt(Math.abs(d));
    ans.push({
        a: (c.x * r + c.y * d) / z,
        b: (c.y * r - c.x * d) / z,
        c: r1,
    });
}

function findTangents(a: Circle, b: Circle): Line[] {
    const ans: Line[] = [];
    // Перебираємо знаки ±r1 та ±r2 (з якого боку проходить пряма)
    for (let i = -1; i <= 1; i += 2)
        for (let j = -1; j <= 1; j += 2)
            tangents(sub(b, a), a.r * i, b.r * j, ans);
    // Зворотний зсув: повертаємо центр першого кола з початку координат
    for (const l of ans)
        l.c -= l.a * a.x + l.b * a.y;
    return ans;
}

Go

func tangents(c Pt, r1, r2 float64, ans *[]Line) {
    r := r2 - r1
    z := sqr(c.X) + sqr(c.Y)
    d := z - sqr(r)
    if d < -EPS {
        return
    }
    d = math.Sqrt(math.Abs(d))
    *ans = append(*ans, Line{
        A: (c.X*r + c.Y*d) / z,
        B: (c.Y*r - c.X*d) / z,
        C: r1,
    })
}

func findTangents(a, b Circle) []Line {
    var ans []Line
    // Перебираємо знаки ±r1 та ±r2 (з якого боку проходить пряма)
    for i := -1; i <= 1; i += 2 {
        for j := -1; j <= 1; j += 2 {
            tangents(Pt{b.X - a.X, b.Y - a.Y}, a.R*float64(i), b.R*float64(j), &ans)
        }
    }
    // Зворотний зсув: повертаємо центр першого кола з початку координат
    for k := range ans {
        ans[k].C -= ans[k].A*a.X + ans[k].B*a.Y
    }
    return ans
}

Задачі

TIMUS 1163 Chapaev