Алгоритмическая задача

Дано: есть матрица высот, то есть (грубо говоря) равномерная сетка точек на поверхности Земли, для каждой из которых приведена ее высота. Точек много, порядка миллиарда. Нужно построить TIN-модель (triangular irregular network), то есть сетку из треугольников с «абы где» расположенными вершинами, «похожую» на исходную матрицу высот. Треугольников должно быть много меньше — несколько десятков тысяч это максимум.

Стандартные алгоритмы построения TIN построены по общему принципу — для каждой точки матрицы высот (всего их N) вычисляется какая-то численная характеристика, связанная с локальными свойствами этой точки, в качестве вершин треугольников выбираются точки с максимальным значением этой характеристики (а их число обозначим M), а затем для построенных вершин строится триангуляция Делоне. Так вот, поговорим про выбор вершин этих треугольников.

В руководстве, откуда я взял это описание, предлагается сначала вычислить характеристику, а уж потом выбрать точки, где она принимает наибольшие значения. Если прикинуть расход памяти для такого действия, то становится очень грустно (не забываем о том, что у нас уже есть загруженная матрица высот). Хочется сделать это несколько более эффективно.

Легко понять, что функция f, выбирающая из последовательности x0, …, xn M максимальных элементов, является индуктивной, то есть f(x0, …, xn) зависит только от f(x0, …, xn-1) и xn. Можно записать это следующим образом: f(x0, …, xn) = F(f(x0, …, xn-1), xn). Индуктивная функция от последовательности длины N вычисляется за N шагов очень простым циклом, надо знать лишь значение f на пустой или одноэлементной последовательности и функцию F.

Очень просто написать функцию, находящую M «наибольших» элементов в N-элементной последовательности за O(N*M) с расходом памяти O(M) — выбранные элементы хранятся в массиве, каждый новый элемент сравнивается со всеми выбранными до этого. Если включить мозг, то можно улучшить оценку до O(N*log M) с той же оценкой расхода памяти (храним выбранные M элементов в самобалансирующемся дереве). А можно ли сделать это еще быстрее, возможно, отказавшись от подхода с индуктивными функциями? Я ответа пока что не знаю.

Алгоритмическая задача: 5 комментариев

  1. Мне вообще не нравится этот твой общий принцип.
    Предположим у тебя правая половина матрицы — параболоид, левая — симметричный параболоид + высокочастотный белый шум с большой амплитудой (возможно сглаженный).
    Насколько я понимаю, твой подход приведёт к тому что все 10k точек будут выбраны из той половины где белый шум: там дохера точек с высокой локальной кривизной, дохера лочек перегиба, в разные стороны скачут все производные, и т.п — какую бы ты метрику ни выбирал, все точки пойдут налево.

    Вот что я придумал за 5 минут:
    Я бы начал с регулярной сетки (например R*R=50×50=2500) — это гарантирует то что макро-вещи будут хоть как-то отражены. Потом например
    1. Сдвигаем каждую точку но не дальше чем на ( sqrt(N) / R ) * S где S например 0.5 или 0.8 чтобы точка располагалась в наиболее экстремальном месте.
    2. Осталось расставить остальные 10k точек. Для этого придумываем как для каждого треугольника оценить отклонение исходных данных от плоскости треугольника, подешевле и приблизительно (думаю честное решение будет довольно дешовым в терминах O(n) и довольно дорогим в миллисекундах/тактах), сравниваем треугольники между собой, досыпаем оставшиеся точки в херовые треугольники.
    3. Долго ебёмся оптимизируя произвидительность: 10E9 точек это очень дохера для таких алгоритмов. Даже для твоего более простого подхода когда все точки идут налево очень дохера.
    4. Долго ебёмся подбирая разные коэффициенты (сколько точек ставить регулярно сколько добавлять потом, насколько далеко можно двигать регулярные точки, чо за характеристику использовать) и тестируя на наборе реальных данных (т.е. в зависимости от требования вычисляя какую-нить статистику нисколько хорошо треугольная сетка сапроксимировала исходный dataset: может надо минимизировать отклонение по высоте, может абсолютное, может разницу объёма, может ещё что)
    5. ???
    6. Profit, бухаем.

    P.S. А чо это за точки вообще? 3D сканер, эхолот с корабля в дно, радар с самолёта/спутника в землю, ещё чо?

    P.P.S. Вариант 2: разбить точки на квадратики например 100×100, для каждого посчитать discrete cosine transform (лучше на GPU — прекрасно для этого подходит), из этих коэффициентов вывести один скаляр на каждый квадратик «шероховатость» т.е. «насколько много тут высоких частот», дальше поделить 10k точек по квадратикам чтобы в каждом квадратике была хотя бы 1 точка, а количество добавочных точек было пропорционально шероховатости.

    1. Идея эта не моя, используется вообще довольно активно, в частности, в ArcGIS, есть там такая функция.

      http://mapcontext.com/autocarto/proceedings/auto-carto-8/pdf/systematic-selection-of-very-important-points(VIP)-from-digital-terrain-models.pdf — вот ссылочка с описанием алгоритма.

      Природа данных — неизвестна, с равным успехом это могут быть как ДЗЗ со спутника, так и сгенерированная по горизонталям топокарты интерполированная картинка.

      Другие подходы тоже посмотрю, но повторюсь, что N непристойно большое.

      1. Откуда у двух чуваков в 1987 году порядка одного гигабайта памяти?

        Природа данных известна: это поверхность планеты Земля. Для алгоритма важно, что на ней нет придуманных мною штук, и присутствуют все компоненты частоты.
        Кстати алгоритм крупно лажает если ему дать поверхность дюны, или местность с замёрзшим озером: будут выбраны точки примерно на границах озера (посередине не будут), но там где берега гладкие они будут сильно выше уровня воды. А если данные восстановлены из изолиний всё будет охуенно для замёрзших океанов (берег=изолиния 0) и лажать для замёрзших озёр и рек.

        Соответственно для искусственных конструкций (поверхность машины яхты самолёта здания) алгоритм непригоден совсем.

        У моего варианта в P.P.S. который вроде не страдает от изложенных недостатков (вероятно страдает от других) сложность по времени { памяти } где-то между O(N + M) { N + M } если размер квадратика константа и O(N*log(N) + M) { M } если кол-во квадратиков константа. И прекрасно переносится на GPU. Можно дальше улучшать (например натравить описанный в статье алгоритм на каждый квадратик по-отдельности, для каждого выбрать и отсортировать 20% точек).

        1. > Откуда у двух чуваков в 1987 году порядка одного гигабайта памяти?

          Заметь, я нигде не говорил о том, что мне понадобится гигабайт памяти. Весь этот миллиард точек может быть записан хоть на стриммере и считываться последовательно. Другой вопрос в том, что вычисление «важности» точки довольно напряжное «в миллисекундах/тактах» и занимает приличное время даже на «недохлой» машине четырехлетней давности.

          > Кстати алгоритм крупно лажает если ему дать поверхность дюны, или местность с замёрзшим озером: будут выбраны точки примерно на границах озера (посередине не будут), но там где берега гладкие они будут сильно выше уровня воды. А если данные восстановлены из изолиний всё будет охуенно для замёрзших океанов (берег=изолиния 0) и лажать для замёрзших озёр и рек.

          Замерзшие моря-океаны, озера и реки восстанавливаются довольно неплохо. Точки на их границе (с нулевой высотой) выбираются часто, так как с одной стороны от них обычно находится склон, а с другой — ровная поверхность. Плюс, довольно часто выбираются точки с «перегиба» склона чуть выше, по-моему, выглядит неплохо.

  2. >можно ли сделать это еще быстрее
    Можно.
    За N+M времени и N памяти.
    Если можно портить массив с исходными данными т.е. они больше не нужны, то за N+M времени с расходом памяти sqrt(N).

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *