Билинейная интерполяция по целочисленным координатам в триангуляции Делоне

У меня есть плоская триангуляция Делоне, состоящая примерно из 1 миллиона треугольников. Каждая вершина помечена несколькими скалярными метриками [1], и я хотел бы увидеть быструю и простую интерполяцию каждой из этих метрик на одной и той же регулярной сетке. Для справки, объединение моих треугольников покрывает около 10 миллионов ячеек сетки, имеющих (целочисленные) координаты. [2]

Когда я говорю просто, я имею в виду просто. Билинейка это хорошо! Насколько я понимаю, это (а) в основном то, чем графические процессоры зарабатывают на жизнь, и (б) вероятно, является предметом бесчисленных домашних заданий. Я сам являюсь государственным исследователем в области общественного здравоохранения, поэтому для меня это не домашнее задание. :-)

В моей медленной, но правильной эталонной реализации я могу вычислить следующее примерно за 10 минут:

Для каждого треугольника T:

  1. Набор G всех (целочисленных) декартовых координат в пределах ограничивающей рамки T;
  2. Барицентрические координаты (u, v, w) для каждого (x, y) в G;
  3. Отказ от (u, v, w), которые не все положительны, то есть внутри T;
  4. Взвешенная сумма (uz_1 + vz_2 + w*z_3) для каждой оставшейся координаты в T, где z_1, z_2 и z_3 — для данной метрики [1] скалярные значения в вершины Т.

Мне действительно нужно, чтобы шаги 1-3 были быстрыми; шаг 4 тривиален, но это моя конечная цель. В идеале решение должно принимать одну из следующих форм:

  • Соответствующая лицензия (GPL в порядке) библиотека с очень простым API; или же
  • Объяснение настолько ясное, что становится очевидным, как программист среднего уровня может написать код на Fortran, R, Python или C.

Классической постановкой этой задачи является работа по моделированию местности «TIN to DEM». Но кажется, что в наши дни чаще требуется обратное (?)

Некоторая базовая очистка, например удаление дубликатов, когда точка попадает точно на ребро или вершину, общую для 2+ треугольников, тоже подойдет.

Заранее большое спасибо за ваше время и внимание. Я уберу форматирование и отредактирую в соответствии с предложениями, как только сойду с поезда!

Сноски:

[1] Высота над уровнем моря, температура и влажность. [2] Целые числа в том смысле, что они расположены на расстоянии 20x20 м друг от друга в сетке UTM. Так что просто масштабируйте на 20.


person dholstius    schedule 12.05.2016    source источник
comment
Быстрые вопросы: какая точность вам нужна для вывода и, если она не чрезмерна, можете ли вы просто использовать, скажем, OpenGL/OpenGL-ES на графическом процессоре? Как вы сказали, это именно то, что делают графические процессоры (и делают это чрезвычайно быстро).   -  person Simon F    schedule 12.05.2016
comment
Мне нужна точность около 1% (возможно, 0,1%) от полной шкалы каждой метрики. Так, например, если метрикой z были температуры, и они колебались от 0 до 100 ° C по всей области, то мне нужна только интерполяция z с точностью около 1,0 или, может быть, 0,1 ° C. Многие треугольники будут плоскими, имеющими вершины. все они очень близки к одному и тому же z. Я попытался сгенерировать триангуляцию таким образом, чтобы градиент был достаточно гладким, то есть там, где я ожидаю, что градиент будет крутым, я создал более мелкую сетку. Я мог бы построить гистограмму диапазона z по размеру треугольника ...?   -  person dholstius    schedule 12.05.2016
comment
Это выглядит многообещающе, за исключением того, что у меня уже есть триангуляция Делоне, которую я хочу использовать:   -  person dholstius    schedule 12.05.2016
comment
Также многообещающе: github.com/geotrellis.   -  person dholstius    schedule 12.05.2016
comment
Похоже на то, что я хочу: docs.scipy.org/doc/scipy-0.14.0/reference/generated/   -  person dholstius    schedule 12.05.2016
comment
Ооо: r-bloggers.com/   -  person dholstius    schedule 12.05.2016
comment
Десять минут кажутся очень медленными. Я ожидаю, что реализация, скомпилированная на C, будет работать максимум через несколько секунд. Какой язык/среду вы используете?   -  person Yves Daoust    schedule 12.05.2016
comment
Я использую R. Возможно, проблема с вводом-выводом памяти (моя наивная реализация использует около 10 ГБ, чего не должно быть в хорошей реализации, поскольку я храню копии промежуточных шагов). Могут быть накладные расходы от переводчика. И т. д. Я использовал пакет profvis, чтобы устранить некоторые очевидные узкие места, и таким образом я сократил время до 10 минут.   -  person dholstius    schedule 13.05.2016
comment
@dholstius: с большей осторожностью вы, вероятно, сможете остаться ниже 1 ГБ.   -  person Yves Daoust    schedule 13.05.2016


Ответы (1)


Хотя я не вижу ничего, что могло бы объяснить медлительность в вашем описании, вот как бы я с этим справился. Основным ингредиентом действительно является «треугольный сканер».

Начните с сортировки трех вершин по оси Y. Это потребует трех сравнений и всего шести возможных конфигураций. Цикл по целочисленным ординатам от верхнего Y до среднего Y, затем от среднего Y до нижнего Y. Для каждой ординаты пересечения с левой и правой сторонами дают вам интервал.

Цикл по целочисленным абсциссам в этом интервале слева направо. Двойной цикл будет посещать только те узлы сетки, которые принадлежат треугольнику.

Вместо использования барицентрических координат вы можете установить уравнение плоскости, т.е. оценить коэффициенты Z = a X + b Y + c и использовать эту формулу для интерполяции. (Вы даже можете вычислять значения постепенно, т. е. Z(X + 1) = Z(X) + a, но для небольшого количества точек на треугольник я не уверен, что это стоит.)

Легко избежать дублирования точек, полагаясь на простое соглашение: создавать только точки, которые попадают на левую сторону треугольника, а не те, которые попадают на правую сторону (они будут получены треугольником справа).

введите здесь описание изображения

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

Общая рабочая нагрузка будет зависеть от количества треугольников, протяженности домена по оси Y и количества покрытых узлов сетки, учитывая несколько арифметических операций для каждого из этих факторов. Для одного миллиона треугольников и десяти миллионов ячеек сетки время работы менее секунды не является маловероятным.

person Yves Daoust    schedule 12.05.2016
comment
Спасибо. Не могли бы вы немного уточнить, как составить уравнение плоскости, учитывая три значения в вершинах? - person dholstius; 13.05.2016
comment
@dholstius: решить простую систему Zi = a Xi + b Yi + c для i=1, 2, 3. - person Yves Daoust; 13.05.2016
comment
Ах я вижу. Спасибо еще раз. - person dholstius; 13.05.2016