Налево или направо: векторное произведение и поворот трёх точек

Если ты знаешь, что такое векторное произведение на плоскости, что значит его знак и как через него найти площадь треугольника и расстояние до прямой, листай до «Теперь сам».

Как думаешь, где точка

Курьер едет по прямой из точки A (1, 1) в точку B (5, 3). Дом клиента стоит в C (2, 4), склад - в D (6, 2). Какой из них слева по ходу движения, а какой справа? А точка (9, 5)?

Ответь себе до чтения - и заодно придумай, как это посчитать программой. Первое, что приходит в голову, - сравнить наклоны: наклон AB равен 2 / 4 = 0,5, наклон AC равен 3 / 1 = 3. Больше - значит, левее? Проверь это правило на дороге, идущей строго вверх: из (2, 1) в (2, 6). Наклон - деление на ноль. А на дороге, идущей справа налево, «больше» вдруг начинает значить «правее».

Этот урок даёт правило без этих дыр. Ответ на вопрос: C слева, D справа, (9, 5) лежит на продолжении дороги.

На что это похоже

Водитель на перекрёстке не измеряет углов. Он смотрит, куда крутить руль: налево, направо или ехать прямо. Для этого не нужно знать, что дорога идёт под углом 26,57 градуса к оси абсцисс, - нужно только сравнить два направления: куда едешь сейчас и куда надо.

Векторное произведение и есть такое сравнение двух направлений. Оно не говорит, на сколько градусов поворачивать, зато честно и без вычислительных хитростей говорит, в какую сторону. В геометрических алгоритмах больше и не нужно: выпуклая оболочка, пересечение отрезков, точка внутри многоугольника - все они спрашивают только сторону.

Механизм: два умножения и одно вычитание

Точка на плоскости - пара чисел (x, y). Вектор - сдвиг от одной точки к другой: из A в B это B − A = (Bx − Ax, By − Ay). Для A (1, 1) и B (5, 3) вектор u = (4, 2): четыре шага вправо, два вверх.

Векторное произведение двух векторов на плоскости - число

u × v  =  ux · vy  −  uy · vx

Его смысл - ориентированная площадь параллелограмма, натянутого на u и v. Модуль - обычная площадь. Знак - в какую сторону поворачивать от u, чтобы прийти к v кратчайшим путём: плюс - против часовой стрелки (налево), минус - по часовой (направо), ноль - векторы лежат на одной прямой.

Откуда здесь именно эта формула, а не скалярное произведение или угол. Повернём u на 90 градусов налево: получится (−uy, ux). Векторное произведение - это скалярное произведение v с этим повёрнутым вектором: −uy · vx + ux · vy. Скалярное произведение p · q = px · qx + py · qy равно произведению длин на косинус угла между векторами, поэтому оно положительно, когда угол острый, - векторы смотрят в одну сторону. Значит, u × v > 0 ровно тогда, когда v смотрит туда же, куда «левая рука» от u.

Для трёх точек берут оба вектора из A:

поворот(A, B, C)  =  (B − A) × (C − A)
                  =  (Bx − Ax)(Cy − Ay)  −  (By − Ay)(Cx − Ax)

Правило

Поворот трёх точек (B − A) × (C − A) положителен, если C слева от направленной прямой A → B (обход A, B, C против часовой стрелки), отрицателен, если справа, и равен нулю, если три точки на одной прямой. Модуль - удвоенная площадь треугольника ABC. В формуле только умножения и вычитания, поэтому на целых координатах она точна.

Разбор: курьер, дом и склад

def povorot(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def skazat(z):
    if z > 0:
        return "налево"
    if z < 0:
        return "направо"
    return "по прямой"

a = [1, 1]
b = [5, 3]
for c in [[2, 4], [6, 2], [9, 5], [3, 2]]:
    z = povorot(a, b, c)
    print(c, z, skazat(z))
print(povorot(b, a, [2, 4]), povorot(a, [2, 4], b))

[2, 4] 10 налево
[6, 2] -6 направо
[9, 5] 0 по прямой
[3, 2] 0 по прямой
-10 -10

Дом (2, 4): 4 · 3 − 2 · 1 = 10, слева. Склад (6, 2): 4 · 1 − 2 · 5 = −6, справа. Точки (9, 5) и (3, 2) дают ноль: обе на прямой AB. Ноль ничего не говорит о том, лежит ли точка между A и B, - это второй вопрос, и следующий разбор его решает.

Последняя строка проверяет, что знак зависит от порядка. Если ехать из B в A, дом оказывается справа: −10. Если поменять местами вторую и третью точку, тоже −10. Отсюда практическое правило: поворот меняет знак при перестановке любых двух точек. Поэтому в каждой функции, которая им пользуется, порядок аргументов - часть смысла, а не деталь.

Вертикальная дорога из (2, 1) в (2, 6), на которой ломалось сравнение наклонов: вектор (0, 5), для точки (1, 3) поворот 0 · 2 − 5 · (−1) = 5 - слева, как и должно быть. Деления нет, особого случая тоже.

Разбор: площадь, расстояние и «лежит ли на отрезке»

Модуль поворота - удвоенная площадь треугольника. Площадь треугольника - это половина основания на высоту, поэтому высота, то есть расстояние от C до прямой AB, равна |поворот| / |AB|. Здесь деление и корень появляются, но только когда нужно само расстояние, а не сторона.

Лежит ли C на ОТРЕЗКЕ AB: поворот равен нулю, и обе координаты C между координатами концов. Второе условие без первого не работает: точка (3, 1) внутри прямоугольника между A и B, но не на отрезке.

import math

def povorot(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def na_otrezke(a, b, c):
    if povorot(a, b, c) != 0:
        return False
    v_x = min(a[0], b[0]) <= c[0] <= max(a[0], b[0])
    v_y = min(a[1], b[1]) <= c[1] <= max(a[1], b[1])
    return v_x and v_y

a = [1, 1]
b = [5, 3]
dlina = math.sqrt((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2)
for c in [[2, 6], [6, 2], [3, 2], [9, 5]]:
    z = povorot(a, b, c)
    print(c, abs(z) / 2, round(abs(z) / dlina, 4), na_otrezke(a, b, c))

[2, 6] 9.0 4.0249 False
[6, 2] 3.0 1.3416 False
[3, 2] 0.0 0.0 True
[9, 5] 0.0 0.0 False

Столбцы: точка, площадь треугольника с A и B, расстояние до прямой, лежит ли на отрезке. (3, 2) и (9, 5) обе на прямой, но на отрезке только первая: у второй x = 9 больше 5. Корень √20 ≈ 4,47 здесь неизбежен - длина AB иррациональна, - но в решение о стороне он не входит.

Что ломается без этого

Наклоны вместо поворота. Деление на ноль на вертикальной прямой, смена смысла «больше» при движении справа налево и дробные числа там, где вход был целым. Каждая из трёх бед порождает свою ветку if, и одна из веток обязательно окажется неверной.

Углы через atan2. Угол по целым координатам обычно иррационален и хранится в float приближённо, и сравнение двух почти равных углов даёт случайный ответ. Поворот на целых координатах даёт ровно 0, когда точки на прямой, и ровно правильный знак, когда нет.

Экранные координаты. На экране, в SVG и в canvas ось y направлена вниз. Та же формула там переворачивается: положительный поворот на экране - это поворот по часовой стрелке. Код, перенесённый из учебника в отрисовку без этой поправки, рисует оболочку «наизнанку».

Зачем это в работе

Навигатор решает, сказать ли «поверните налево»: поворот текущего отрезка маршрута и следующего. Сервис доставки проверяет, по какую сторону от реки или кольцевой дороги стоит адрес. Видеокарта отбрасывает грани, повёрнутые от камеры, по знаку того же числа, посчитанного для треугольника на экране. Библиотеки геометрии (CGAL, GEOS под Shapely и PostGIS) строят проверки сторон и пересечений поверх функции, которая называется orient2d или orientation, - это она (проверка «в окружности» стоит на её родственнице incircle).

Когда ты видишь геометрический алгоритм на тысячу строк, полезно найти в нём это место: почти всегда решения принимаются в пяти-шести вызовах поворота, а остальное - перебор кандидатов.

Теперь сам

1. A (0, 0), B (4, 1). Слева или справа точки (1, 3), (3, −2), (8, 2)?
2. Площадь треугольника (0, 0), (6, 0), (2, 5) - через поворот.
3. Чему равен поворот(A, A, C) и что это значит для алгоритма, получившего две одинаковые точки?