Площадь многоугольника

Сколько квадратных километров в Медеуском районе Алматы? Акимат района пишет: 253,4. Когда в 2023 году утверждали новый генеральный план, площадь всего города назвали такую: 70 348 гектаров, то есть 703,48 км². А если сложить площади восьми районов, которые приводит Википедия со ссылками на сайты акиматов, выйдет 684,5 км²: девятнадцать квадратных километров куда-то делись. Проверить можно самим. В OpenStreetMap у каждого района есть граница — замкнутая ломаная из нескольких сотен точек, и площадь по ним считает формула, которая обходит границу и складывает треугольники. Выберите район и посмотрите, как она это делает.

Синие треугольники прибавляются, оранжевые вычитаются, и всё, что лежит вне района, сокращается — куда бы вы ни поставили точку O. Для семи районов из восьми ответ расходится с опубликованным не больше чем на 1,1%: граница в OpenStreetMap нарисована не идеально, а год, к которому относятся опубликованные цифры, не указан. Медеуский район выбивается: 273,22 км² против 253,4 — на 19,8 км² больше. При этом сумма всех восьми районов по OpenStreetMap — 702,91 км², на 0,08% меньше площади города по генплану. Так «потерянные» девятнадцать квадратных километров нашлись: они в Медеуском районе. Граница в OpenStreetMap и площадь города их учитывают, а опубликованная площадь района, видимо, нет. Почему цифры разошлись, по открытым источникам выяснить не удалось — но заметить расхождение помогла одна строка кода.

Идея

Многоугольник задан вершинами P₀, P₁, …, P_n−1 в порядке обхода. Возьмём любую точку O и соединим её со всеми вершинами: получится веер из n треугольников OP_iP_i+1, по одному на сторону. Площадь каждого даёт векторное произведение: половина cross(O, P_i, P_i+1), со знаком плюс, если из точки O сторона видна идущей против часовой стрелки, и минус — если по часовой.

Площадь многоугольника — сумма площадей со знаком всех треугольников веера, из какой бы точки O он ни был построен. Если O — начало координат, треугольник на стороне P_iP_i+1 даёт ½(x_iy_i+1 − x_i+1y_i), и весь алгоритм — одна сумма по сторонам.

Это формула шнурков. Без модуля сумма — площадь со знаком: плюс, если вершины перечислены против часовой стрелки, минус — если по часовой. Формула верна для простого многоугольника — такого, у которого стороны не пересекаются и не касаются, кроме соседних в общей вершине. Выпуклость не нужна: район может быть каким угодно извилистым.

Почему шнурки

Выпишите координаты вершин столбиком и повторите первую вершину внизу. Каждое слагаемое — крест между соседними строками: x_i по диагонали вниз к y_i+1 — со знаком плюс, y_i вниз к x_i+1 — со знаком минус. Перекрещённые линии похожи на шнуровку ботинка, отсюда английское название shoelace formula. Тащите вершины, обходите многоугольник в другую сторону, сделайте из него бантик.

Стрелка из пяти вершин — тот же пример, что в листинге ниже. Слагаемые −5, 30, 0, 20 и −5 дают в сумме 40, площадь — 20 клеток. Ноль получается у стороны, которая лежит на прямой, проходящей через начало координат: её треугольник вырожден. Отрицательные слагаемые не ошибка — это треугольники, которые нужно вычесть. Если обойти стрелку в другую сторону, сумма станет −40: формула замечает направление обхода.

У той же суммы есть вторая запись — формула трапеций. Опустите из концов каждой стороны перпендикуляры на ось x: получится трапеция с площадью ½(x_i − x_i+1)(y_i + y_i+1). При обходе против часовой стрелки стороны верхнего края фигуры идут справа налево и дают трапеции со знаком плюс, стороны нижнего края — слева направо и со знаком минус. Разность оставляет саму фигуру. Раскройте скобки: произведения x_iy_i при суммировании по кругу сокращаются, и остаётся формула шнурков.

Почему это площадь

Сначала убедимся, что ответ не зависит от выбора точки O. Сдвинем O в точку v. Каждое слагаемое изменится: cross(P_i − v, P_i+1 − v) = cross(P_i, P_i+1) − cross(P_i, v) + cross(P_i+1, v). Но при суммировании по кругу добавки сокращаются: каждая вершина входит в одну добавку со знаком плюс, а в соседнюю — с минусом. Сумма та же. Поэтому в сцене можно таскать O куда угодно, и итог не меняется.

Сцена показывает ту же теорему с другой стороны. Возьмите любую точку плоскости и посчитайте, сколько синих и сколько оранжевых треугольников веера её накрывают. Внутри района одних окажется на один больше, чем других, снаружи — поровну. Граница многоугольника — это сумма границ треугольников веера: спицы от O к вершинам проходятся дважды, туда и обратно, и сокращаются. Поэтому разность «синие минус оранжевые» для точки равна тому, сколько раз граница обходит её против часовой стрелки: внутри простого многоугольника это плюс или минус единица, смотря по направлению обхода, снаружи — ноль. Это дискретный вариант формулы Грина из анализа, где площадь области выражается интегралом по её границе.

Знак — это направление обхода

Без модуля формула отвечает на второй вопрос: в какую сторону перечислены вершины. Плюс — против часовой стрелки, минус — по часовой. В OpenStreetMap направление линий границы ничего не значит, так прямо и сказано в правилах для мультиполигонов, и у районов Алматы оно разное: у Алмалинского, Жетысуского и Наурызбайского границы записаны против часовой стрелки, у остальных — по часовой.

Форматы, где направление важно, договорились о нём по-разному. В GeoJSON (RFC 7946, 2016 год) внешнее кольцо многоугольника обходится против часовой стрелки, а дыры — по часовой. В шейп-файлах ESRI (описание формата 1998 года) наоборот: внутренность должна оставаться справа от идущего по кольцу, так что внешние кольца идут по часовой стрелке, а дыры — против. Договорённость нужна ради дыр: если кольца направлены по правилу, площадь фигуры с дырами — просто сумма площадей колец со знаком, без единого модуля, и знак каждого кольца заодно говорит, дыра оно или нет.

А если не направлены, сумма со знаком врёт. Медеуский район в OpenStreetMap — два кольца: основная часть и микрорайон Алатау в восьми с половиной километрах к северо-востоку от неё. Основная часть записана по часовой стрелке и даёт −267,38 км², Алатау — против часовой и даёт +5,85 км². Сумма со знаком — −261,53 км², а настоящая площадь — сумма модулей, 273,22 км². Выберите в первой сцене Медеуский район и дождитесь конца обхода. Когда направление колец не гарантировано, модуль берут у каждого внешнего кольца отдельно, а дыры вычитают по их роли в данных, а не по знаку. Об этом первая задача.

Код

Формула — одна строка. Листинг считает её в целых числах и рядом показывает вариант с трапециями, центр тяжести и то, что бывает с самопересекающимся бантиком. Последние строки — ловушка с плавающей точкой, о ней ниже.

Центр тяжести считается той же суммой: каждый треугольник веера весит своей площадью со знаком, а его центр — среднее трёх вершин. Так первая сцена ставит подписи районов. У бантика четыре вершины и пересекающиеся стороны. Формула даёт ноль, хотя бантик закрывает восемь клеток: его половины обходятся в разные стороны и вычитаются. Сумма шнурков честно складывает площадь с учётом того, сколько раз и в какую сторону граница обходит каждую точку, — но для самопересекающейся границы это не та площадь, которую хотелось. Поэтому данные, где такое возможно, сначала проверяют на пересечения сторон.

А это первая сцена в цифрах: восемь районов из OpenStreetMap, их площади рядом с опубликованными и Медеуский район по кольцам.

Ловушка: градусы — не метры

Координаты в OpenStreetMap — широта и долгота в градусах. Подставьте их в формулу как есть — получится площадь в «квадратных градусах»: у Турксибского района 0,008378. Одним множителем её в километры не перевести. Градус широты везде около 111 км, а градус долготы на широте Алматы — около 81 км и к северу короче: меридианы сходятся к полюсу. Последний блок листинга считает площадь одного и того же Турксибского района на трёх картах в метрах:

Turksib, one polygon, four maps:
  degrees as they are:          0.008378 square degrees
  the road map's grid:          75.12 km2
  one scale at 43.2425 N:       75.59 km2
  cos(latitude) at each vertex: 75.47 km2

Сетка карты дорог. На ней стоят сцены каталога с улицами Алматы, и для градуса широты она берёт 110 574 м — длину градуса у экватора. У Алматы градус широты на полпроцента длиннее, 111 097 м. Для рисования это неважно, а площадь выходит меньше: 75,12 км².

Один масштаб на широте центра карты, 43,2425°. Это честные длины градусов, но взятые в одной точке. Турксибский район лежит севернее, где градус долготы короче, и единый масштаб завышает площадь на 0,16%: 75,59 км².

Свой масштаб у каждой вершины: метры на восток — это разница долгот, умноженная на радиус и на косинус широты этой вершины. Такая карта называется синусоидальной проекцией, и площади на шаре она сохраняет. Маленькая клетка размером dφ по широте и dλ по долготе имеет на шаре площадь R² cos φ dφ dλ. На карте она сдвигается вбок и становится параллелограммом, но высота у него R dφ, а основание R cos φ dλ — площадь та же. Это 75,47 км². Для проверки мы пересчитали все восемь районов по эллипсоиду WGS 84 библиотекой GeographicLib: расхождение не больше 0,003%.

Для города этого хватает. Для страны или океана — нет: стороны многоугольника на Земле — дуги, и считать надо на сфере или эллипсоиде. Но и там работает та же идея. GeographicLib Чарльза Карни складывает по каждой стороне площадь между ней и экватором и получает площадь многоугольника — трапеции землемера, только на эллипсоиде.

Узлы сетки: формула Пика

Если все вершины лежат в узлах целочисленной сетки, у площади появляется второй смысл. Формула Пика: S = I + B/2 − 1, где I — узлы строго внутри, B — узлы на границе вместе с вершинами. Сцена со шнурками показывает оба числа: у стрелки 18 узлов внутри и 6 на границе, 18 + 3 − 1 = 20 — та же площадь, что дала сумма шнурков. Доказательство разрезанием на треугольники, с картинками, — в главе «Царицы наук» о площади и объёме.

Программисту формула Пика полезна в обратную сторону. Площадь известна точно из шнурков. Узлы на границе считаются без перебора: на стороне с концами в узлах их gcd(|Δx|, |Δy|) + 1, а наибольший общий делитель даёт алгоритм Евклида. Из формулы выходит число узлов внутри — без единой проверки «лежит ли точка внутри». В Ауэзовском районе, переведённом в целые метры, так получается 23 481 364 узла метровой сетки; это третий тест задачи о яблонях ниже.

Георг Пик опубликовал формулу в 1899 году. Широко известной её сделала книга Гуго Штейнгауза «Математический калейдоскоп» — в английском издании 1950 года. Сам Пик этого признания не застал: в июле 1942 года его, восьмидесятидвухлетнего, депортировали в концлагерь Терезиенштадт, где он через две недели погиб.

Сколько это стоит

Один проход по вершинам, одно умножение и вычитание на сторону: O(n) времени и O(1) памяти сверх самих вершин. Быстрее нельзя: от каждой вершины зависит ответ, и её надо хотя бы прочитать. Восемь районов Алматы — это 7158 вершин, то есть 7158 умножений-вычитаний.

Интересное в этой сложности — точность. Для целых координат ответ точен: удвоенная площадь — целое число. В Go стоит подумать о переполнении. При координатах до 10⁹ удвоенная площадь не больше 8·10¹⁸ и помещается в int64, а промежуточные суммы могут выйти за предел. Но знаковые целые в Go при переполнении считаются по модулю 2⁶⁴, это закреплено в спецификации языка, и итог, который сам помещается в int64, выходит верным. Python о переполнении не думает вовсе.

С плавающей точкой хуже, и подвох — в расстоянии до начала координат. Последние строки первого листинга сдвигают стрелку площадью 20 на 10⁷, 10⁸ и 10⁹ по обеим осям и считают её в числах double. При сдвиге на 10⁸ удвоенная площадь выходит 38 вместо 40, при сдвиге на 10⁹ — ноль. Произведения координат около 10¹⁸, а double хранит целые точно только до 2⁵³ ≈ 9·10¹⁵, и разность двух огромных почти равных произведений теряет все значащие цифры. Координаты OpenStreetMap в единицах 10⁻⁷ градуса — как раз числа такого порядка. Лекарство простое: строить веер из первой вершины, а не из начала координат. Вычтите P₀ из всех точек, числа станут маленькими, и листинг при любом из трёх сдвигов печатает 40.

Где работает

• Землемерие и кадастр. Площадь участка по координатам межевых знаков. Барт Брейден в 1986 году писал, что формула хорошо известна землемерам, но её почти нет в учебниках математического анализа, — и предлагал это исправить.
• Геоданные. Площади районов, озёр, лесов, зданий по их границам: в проекции, сохраняющей площади, или на эллипсоиде, как в GeographicLib.
• Компьютерная графика. Спецификация OpenGL решает, повёрнут ли треугольник к зрителю лицом или спиной, по знаку его площади в оконных координатах и записывает эту площадь той же суммой шнурков. Если включено отсечение граней, треугольники, повёрнутые спиной, отбрасываются и вовсе не закрашиваются.
• Форматы геоданных. Чтобы записать многоугольник в GeoJSON или шейп-файл по правилам формата, направление каждого кольца проверяют по знаку площади и при необходимости переворачивают кольцо.
• Центр тяжести. Подпись района на карте, центр участка, точка опоры плоской детали — та же сумма с весами.
• Олимпиады. Площадь по координатам, направление обхода и формула Пика — частые гости задач по геометрии, особенно там, где нужен точный ответ в целых числах.

История

Формулу описал в 1769 году профессор Альбрехт Людвиг Фридрих Мейстер, который среди прочего преподавал землемерие. Её называют также формулой Гаусса и формулой землемера, а английское прозвище shoelace — за перекрещённые произведения. В школьной геометрии площадь многоугольника ищут разрезанием на треугольники и прямоугольники. Формула по координатам, где разрезания не видно вовсе, делает то же самое: треугольники веера и есть разрезание, только со знаками.

Задачи

В первой задаче участки состоят из нескольких колец с любым направлением обхода, как Медеуский район, а координаты доходят до миллиарда, где double уже ошибается. Во второй — яблони в узлах сетки: формула Пика в обратную сторону, с наибольшими общими делителями на сторонах.