Ответ 1
Я думаю, что следующее может быть полезным приближением, которое масштабируется линейно, а не квадратично с числом точек, и его довольно легко реализовать:
- вычислить центр масс M точек
- найдите точку P 0, которая имеет максимальное расстояние до M
- найдите точку P 1, которая имеет максимальное расстояние до P 0
- приблизительный максимальный диаметр с расстоянием между P 0 и P 1
Это можно обобщить, повторив шаг 3 N раз, и расстояние между P N-1 и P N
Шаг 1 может быть эффективно использован для приближения М к средним значениям долгот и широт, что хорошо, когда расстояния "малы", а полюса достаточно далеко. Другие этапы могут быть выполнены с использованием точной формулы расстояния, но они намного быстрее, если координаты точек можно аппроксимировать как лежащие на плоскости. Как только "далекая пара" (надеюсь, пара с максимальным расстоянием) была найдена, ее расстояние можно пересчитать с помощью точной формулы.
Примером аппроксимации может быть следующее: если φ (M) и λ (M) - широта и долгота центра масс, рассчитанная как Σφ (P)/n и Σλ (P)/n,
- x (P) = (λ (P) - λ (M) + C) cos (φ (P))
- y (P) = φ (P) - φ (M) [это только для ясности, оно также может быть просто y (P) = φ (P)]
где C обычно 0, но может быть ± 360 °, если набор точек пересекает линию λ = ± 180 °. Чтобы найти максимальное расстояние, вам просто нужно найти
- max ((x (P N)) - x (P N-1)) 2 + (y (P N) - y (P N-1)) 2)
(вам не нужен квадратный корень, потому что он монотонен)
Такое же преобразование координат можно было бы использовать для повторения шага 1 (в новой системе координат), чтобы иметь лучшую начальную точку. Я подозреваю, что если выполняются некоторые условия, вышеуказанные шаги (без повторения шага 3) всегда приводят к "истинной далекой паре" (моя терминология). Если бы я только знал, какие условия...
EDIT:
Я ненавижу строить решения других, но кому-то придется.
Сохраняя вышеуказанные 4 шага, с необязательным (но, вероятно, полезным, в зависимости от типичного распределения точек) повторением шага 3, и после решения Spacedman, выполнение вычислений в 3D преодолевает ограничения близости и расстояния от полюсов:
- x (P) = sin (φ (P))
- y (P) = cos (φ (P)) sin (λ (P))
- z (P) = cos (φ (P)) cos (λ (P))
(единственное приближение состоит в том, что это справедливо только для идеальной сферы)
Центр масс определяется как x (M) = Σx (P)/n и т.д. и максимальный, который нужно искать,
- max ((x (P N)) - x (P N-1)) 2 + (y (P N) - y (P N-1)) 2 + (z (P N)) - z (P N-1суб > )) 2)
Итак: сначала вы преобразовываете сферические в декартовы координаты, затем начинаете с центра масс, чтобы найти, по крайней мере, два шага (шаги 2 и 3), самую дальнюю точку из предыдущей точки. Вы можете повторить шаг 3, пока расстояние увеличивается, возможно, с максимальным количеством повторений, но это не приведет вас к локальному максимуму. Исход из центра масс также не очень помогает, если точки распределены по всей Земле.
ИЗМЕНИТЬ 2:
Я достаточно узнал R, чтобы записать ядро алгоритма (хороший язык для анализа данных!)
Для плоского приближения, игнорируя проблему вокруг линии λ = ± 180 °:
# input: lng, lat (vectors)
rad = pi / 180;
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i = which.max((x - mean(x))^2 + (y )^2)
j = which.max((x - x[i] )^2 + (y - y[i])^2)
# output: i, j (indices)
На моем компьютере требуется меньше секунды, чтобы найти индексы i и j для 1000000 точек.
Следующая трехмерная версия немного медленнее, но работает для любого распределения точек (и не необходимо изменить в случае пересечения линии λ = ± 180 °):
# input: lng, lat
rad = pi / 180
x = sin(lat * rad)
f = cos(lat * rad)
y = sin(lng * rad) * f
z = cos(lng * rad) * f
i = which.max((x - mean(x))^2 + (y - mean(y))^2 + (z - mean(z))^2)
j = which.max((x - x[i] )^2 + (y - y[i] )^2 + (z - z[i] )^2)
k = which.max((x - x[j] )^2 + (y - y[j] )^2 + (z - z[j] )^2) # optional
# output: j, k (or i, j)
Вычисление k может быть опущено (т.е. результат может быть задан i и j), в зависимости от данных и требований. С другой стороны, мои эксперименты показали, что вычисление дальнейшего индекса бесполезно.
Следует помнить, что в любом случае расстояние между результирующими точками является оценкой, которая является нижней границей "диаметра" множества, хотя очень часто это будет сам диаметр (как часто это зависит от данные.)
ИЗМЕНИТЬ 3:
К сожалению, относительная ошибка плоского приближения может в крайних случаях достигать 1-1/√3 ≅ 42,3%, что может быть неприемлемым даже в редких случаях. Алгоритм может быть изменен, чтобы получить верхнюю границу приблизительно 20%, которую я получил компасом и прямым фронтом (аналитическое решение громоздко). Измененный алгоритм находит пару точек с локально максимальным расстоянием, а затем повторяет те же шаги, но на этот раз, начиная с середины первой пары, возможно, найдя другую пару:
# input: lng, lat
rad = pi / 180
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i.n_1 = 1 # n_1: n-1
x.n_1 = mean(x)
y.n_1 = 0 # = mean(y)
s.n_1 = 0 # s: square of distance
repeat {
s = (x - x.n_1)^2 + (y - y.n_1)^2
i.n = which.max(s)
x.n = x[i.n]
y.n = y[i.n]
s.n = s[i.n]
if (s.n <= s.n_1) break
i.n_1 = i.n
x.n_1 = x.n
y.n_1 = y.n
s.n_1 = s.n
}
i.m_1 = 1
x.m_1 = (x.n + x.n_1) / 2
y.m_1 = (y.n + y.n_1) / 2
s.m_1 = 0
m_ok = TRUE
repeat {
s = (x - x.m_1)^2 + (y - y.m_1)^2
i.m = which.max(s)
if (i.m == i.n || i.m == i.n_1) { m_ok = FALSE; break }
x.m = x[i.m]
y.m = y[i.m]
s.m = s[i.m]
if (s.m <= s.m_1) break
i.m_1 = i.m
x.m_1 = x.m
y.m_1 = y.m
s.m_1 = s.m
}
if (m_ok && s.m > s.n) {
i = i.m
j = i.m_1
} else {
i = i.n
j = i.n_1
}
# output: i, j
3D-алгоритм может быть изменен аналогичным образом. Возможно (как в 2D, так и в 3D-случае) снова начать с середины второй пары точек (если найдено). Верхняя граница в этом случае "оставлена как упражнение для читателя": -).
Сравнение модифицированного алгоритма с (слишком простым) алгоритмом показало, что для нормального и для квадратного равномерного распределения было почти удвоение времени обработки и уменьшение средней ошибки от 0,6% до 0,03% (порядок). Дальнейший перезапуск из середины приводит к слегка более средней средней ошибке, но почти равной максимальной ошибке.
РЕДАКТИРОВАТЬ 4:
Мне еще предстоит изучить эту статью, но похоже, что 20%, которые я нашел с компасом и прямолинейным, на самом деле 1 -1/√ (5-2√3) ≅ 19.3%
