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

«Ну-ну-ну, что у нас тут? Кажется, мы наткнулись на проблему, требующую решения — как найти расстояние между любыми двумя заданными точками в заданном наборе функций GeoJSON с помощью JavaScript. Но не бойтесь, храбрые рыцари программирования, потому что у нас есть решение, которое сделает эту задачу легкой, как шутка шута. Так что собирайтесь, и давайте отправимся в путешествие математического веселья и славы GeoJSON!»

По словам посредственного заместителя главного инженера, «не великий, не ужасный».

Мотивация статьи — функция, которую мы недавно внедрили на CodeComply.ai. Рассматриваемая функция требует расчета максимального расстояния между любыми двумя точками на заданном плане этажа. Обычно мы работаем с объектами плана этажа, представленными в виде функций GeoJSON, но поскольку визуализировать фактические географические данные намного проще, я попытаюсь измерить самое длинное диагональное расстояние внутри случайной страны, скажем, Польши!

К счастью, в репозитории polska-geojson Piotr Patrzyk есть отличное представление Польши в формате GeoJSON. Давайте использовать это!

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

  <head>
    <link
      rel="stylesheet"
      <link
      rel="stylesheet"
      href="https://unpkg.com/[email protected]/dist/leaflet.css"
      integrity="sha256-kLaT2GOSpHechhsozzB+flnD+zUyjE2LlfWPgU04xyI="
      crossorigin=""
    />

    <script
      src="https://unpkg.com/[email protected]/dist/leaflet.js"
      integrity="sha256-WBkoXOwTeyKclOHuWtc+i2uENFpDZ9YPdf5Hf+D7ewM="
      crossorigin=""
    ></script>
  </head>

  <body>
    <div id="map"></div>
  </body>

Затем мы можем приступить к настройке карты в JavaScript:

import L from "leaflet";

const map = L.map("map").setView([52.069330831388974, 19.480308502370697], 5);

Любознательный читатель может спросить — что это за координаты? Так получилось, что это геометрический центр Польши, который, в свою очередь, оказался в городе, чье название переводится как «пятница», который не является ни серединой рабочей недели, ни нормальной точкой.

Это не впечатляет, и чтобы это стало действительно впечатляющим, нам нужны тайлы карты. Используем OpenStreetMap:

import L from "leaflet";

const map = L.map("map").setView([52.069330831388974, 19.480308502370697], 5);

L.tileLayer("https://tile.openstreetmap.org/{z}/{x}/{y}.png", {
  maxZoom: 19,
  attribution:
    '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>'
}).addTo(map);

Так-то лучше. Далее давайте загрузим данные GeoJSON. Для этого нам нужен слой GeoJSON на нашей карте:

import L from "leaflet";

const map = L.map("map").setView([52.069330831388974, 19.480308502370697], 5);

L.tileLayer("https://tile.openstreetmap.org/{z}/{x}/{y}.png", {
  maxZoom: 19,
  attribution:
    '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>'
}).addTo(map);

const geoJSONLayer = L.geoJSON().addTo(map);

Затем, предполагая, что мы скопировали данные Польши в файл data.json рядом, мы можем добавить его на карту:

import L from "leaflet";
import data from "./data.json";

const map = L.map("map").setView([52.069330831388974, 19.480308502370697], 5);

L.tileLayer("https://tile.openstreetmap.org/{z}/{x}/{y}.png", {
  maxZoom: 19,
  attribution:
    '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>'
}).addTo(map);

const geoJSONLayer = L.geoJSON().addTo(map);

geoJSONLayer.addData(data);

Теперь у нас есть карта с 16 функциями GeoJSON, представляющими польские воеводства. Наши объекты содержат ровно 8192 точки, число, сравнимое с количеством точек, которое мы обычно видим на плане этажа среднего размера. Приступим к нахождению наибольшего диагонального расстояния.

Каким будет самый простой способ сделать это? Как упоминалось ранее, у нас есть набор признаков — многоугольники, каждый из которых представляет собой набор точек. Самым простым вариантом было бы перебрать все точки и для каждой перебрать все остальные, вычислить все расстояния и найти самые длинные. Давайте попробуем это! Во-первых, нам нужно извлечь все точки. Обычно нам нужно извлечь их из свойства geometry.coordinates каждой функции, но поскольку мы используем Turf.js, мы можем использовать помощник, который сделает это за нас:

import { coordAll } from "@turf/meta";
import data from "./data.json";


const points = coordAll(data);

На M1 Max MacBook Pro это занимает около половины миллисекунды для 8192 точек, что незначительно.

Чтобы вычислить наибольшее расстояние, нам сначала нужно вычислить любое расстояние между двумя точками. Здесь необходимо сделать важное примечание: изначально мы рассчитывали расстояния на плане этажа, где на практике кривизна Земли не имеет значения. Делая то же самое на реальных географических объектах, мы должны учитывать это и использовать, например, помощник Turf.js, который использует формулу Хаверсина и работает намного, намного, намного медленнее. Поскольку мы действительно хотим быть быстрыми, мы предполагаем, что Земля плоская (о боже) и используем Теорему Пифагора для расчета расстояния, дешевого и грязного:

const calculatePointsDistance = ([x1, y1], [x2, y2]) => {
  return Math.sqrt(Math.pow(x1 - x2, 2) + Math.pow(y1 - y2, 2));
};

Затем мы можем перебрать точки:

const { maxDistance, p1, p2 } = points.reduce(
  (acc, curr, _, array) => {
    for (let i = 0; i < array.length; i++) {
      const distance = calculatePointsDistance(curr, array[i]);

      if (distance > acc.maxDistance) {
        acc.maxDistance = distance;
        acc.p1 = curr;
        acc.p2 = array[i];
      }
    }

    return acc;
  },
  {
    maxDistance: 0,
    p1: null,
    p2: null
  }
);

В результате вызова reduce мы получаем maxDistance (что само по себе не имеет особого смысла, но является каким-то представлением расстояния) и координирует две точки, которые разделяет максимальное расстояние.

Поскольку вышеприведенная реализация O(n²), и мы запускаем ее на 8192 записях, для запуска требуется около 430 мс, что очень много, если вы хотите сделать это в режиме реального времени в браузере. .

Есть простое улучшение, которое мы можем выполнить — половина итераций, которые мы делаем выше, являются повторениями. Это происходит потому, что мы запускаем проверку декартова произведения массива точек на себя. Это означает, что для каждой пары точек p1, p2 мы делаем такую ​​же проверку и для p2, p1, что не имеет особого смысла. Убираем повторы:

const { maxDistance, p1, p2 } = points.reduce(
  (acc, curr, index, array) => {
    for (let i = index + 1; i < array.length; i++) {
      const distance = calculatePointsDistance(curr, array[i]);

      if (distance > acc.maxDistance) {
        acc.maxDistance = distance;
        acc.p1 = curr;
        acc.p2 = array[i];
      }
    }

    return acc;
  },
  {
    maxDistance: 0,
    p1: null,
    p2: null
  }
);

Это не стало неожиданностью; это занимает примерно половину того, что делала предыдущая реализация, около 220 мс. Что все еще слишком, слишком долго. Но эй, это все еще алгоритм O(n²/2).

Чтобы наше решение работало быстрее, нам действительно нужно подумать о том, что заставляет его работать так долго. Удивительно, но не факт, что проверка, которую мы выполняем, O(n²/2). Это из-за количества баллов! Если вы посмотрите на карту, то увидите, что многие точки на самом деле не нужны — в первую очередь на ум приходят границы между воеводствами. Есть простой способ удалить эти точки — мы можем объединить наши полигоны с помощью помощника Turf.js, который берет два полигона и возвращает их объединение:

import union from "@turf/union";

const mergedPolygon = data.features.reduce((acc, curr) => {
  if (!acc) return curr;

  return union(acc, curr);
}, null);

const mergedPoints = coordAll(mergedPolygon);

Затем, после применения подхода без повторения к новому набору точек (всего 1699 из тех, что остались после слияния полигонов), мы видим, что это примерно на 45% быстрее (это занимает около 120 мс). . Лучше, но все еще нет сигары.

Хорошо, но нужны ли нам все эти точки для расчета максимального расстояния? Рассмотрим банан:

Интуиция подсказывает нам, что некоторые точки, составляющие форму этого 2D-представления банана, на самом деле не имеют значения. На самом деле, может быть, мы могли бы обернуть банан (каламбур) в какую-то упрощенную форму, которая все же позволила бы нам точно рассчитать максимальное диагональное расстояние? Начнем с самых простых вариантов.

Первый заключается в вычислении ограничивающей рамки. Это, однако, будет работать только в случаях, когда максимальное расстояние по диагонали на самом деле является диагональю полученного прямоугольника, так что это слишком специфично для применения. Мы также могли бы попытаться построить наименьший круг, окружающий банан, но это не имеет никакого смысла. К счастью, есть объемная форма, которая нам и нужна: выпуклая оболочка. Почему выпуклый корпус, а не вогнутый? Ну, снова рассмотрим банан. Интуиция подсказывает нам, что точки на внутренней кривой банана не имеют значения. Если это так, то нашему корпусу не нужна вогнутость, которая на самом деле требует больше вычислительной мощности для расчета.

Давайте посчитаем выпуклую оболочку Польши, используя подходящий помощник Turf.js:

import { featureCollection, point } from "@turf/helpers";
import convex from "@turf/convex";

// convex accepts a featureCollection made up of actual GeoJSON points
const collection = featureCollection(points.map((p) => point(p)));
const convexHull = convex(collection);

const convexPoints = coordAll(convexHull);

Теперь все, что нам нужно сделать, это перебрать новый набор точек без повторений. Результат? Что ж, на самом деле это неплохо, так как это занимает 13 % времени метода «объединения полигонов» и 3,7 % исходного подхода грубой силы при 16 мс! Это вполне приемлемо с точки зрения запуска в браузере.

Почему это намного быстрее, чем метод грубой силы? Потому что из исходных 8192 точек выпуклой оболочке нужно только 40, чтобы сохранить охватывающую форму. Почему расчет выпуклой оболочки выполняется быстрее алгоритма O(n²/2)? Ну, потому что на самом деле это ближе к O(n log n), и это в худшем случае! Алгоритму нужно только проверить, находится ли текущая итерируемая точка внутри промежуточного полигона, и если нет — расширить полигон.

Для дальнейшего чтения также стоит упомянуть, что количество вершин различается в зависимости от типа распределения рассматриваемых точек. Например, для распределения Гаусса число очень, очень маленькое — это будет O(√(log n))! Еще один интересный факт заключается в том, что подход с выпуклой оболочкой не самый эффективный, но он простой и достаточно хороший для этого варианта использования.

Хорошо, кажется, что наша работа здесь сделана! Так ли это? Какое реальное расстояние? Как это выглядит на карте? Давайте посмотрим. Во-первых, нам нужно немного изменить карту. Давайте добавим немного цвета слою GeoJSON:

const geoJSONLayer = L.geoJSON(null, {
  style: (feature) => {
    if (feature.geometry.type === "LineString") {
      return {
        color: "red"
      };
    }
  }
}).addTo(map);

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

import { lineString } from "@turf/helpers";
import length from "@turf/length";

// ...

const line = lineString([convexP1, convexP2]);
const len = length(line);

geoJSONLayer.addData(line);

Ну а расстояние (без учета кривизны Земли) 768 с четвертью километра. Что касается графического представления:

РЕДАКТИРОВАТЬ: Основываясь на комментарии на HackerNews, есть еще один способ ускорить расчет: удалить вызов квадратного корня из функции, которая вычисляет расстояние между двумя точками внутри цикла. Дело в том, что на самом деле это не вычисление расстояния; это просто функция, вычисляющая синтетическое значение, которое можно использовать для сравнения пар точек с точки зрения расстояния между ними. Он даже не используется для расчета фактического расстояния. Таким образом, вызов Math.sqrt на самом деле ни для чего не используется, и он удивительно дорог. Его удаление сократило время выполнения примерно до 10 мс, что является огромным преимуществом!

Код для всего вышеперечисленного можно найти в этой песочнице CodeSandbox. Чтобы узнать точное время на своем компьютере, обязательно проверьте панель консоли браузера, а не панель CodeSandbox!

EDIT2: Как было сказано в начале, я использовал евклидов подход для простоты, но также довольно легко получить результаты, учитывающие кривизну Земли. Вместо использования функции для вычисления значения на основе расстояния (для чего я использовал сумму квадратов разностей между координатами) мы можем использовать хелпер расстояния Turf.js, который учитывает кривизну и вычисляет расстояние по формуле гаверсина. Этот подход немного менее эффективен (для той же настройки, что и выше, это занимает около 16 мс), но он дает результаты, которые действительно имеют смысл для сферической карты.

Результаты сильно отличаются от приведенных выше — нижняя правая точка совершенно другая (как показано на карте ниже), а расстояние на сфере больше, чем расстояние марафона! Он длиннее почти на 43 километра — 811 240 км.

Код можно найти в this CodeSandbox.

Я очень надеюсь, что вам понравилась эта статья. Подписывайтесь на меня, чтобы не пропустить другие статьи о геометрии, географии, компьютерном зрении, оптимизации и обо всем этом (и многом другом) в JavaScript!

Я наставляю разработчиков программного обеспечения. Напишите мне на MentorCruise для долгосрочного наставничества или на CodeMentor для индивидуальных сессий.