Задача интерполяции, в общем случае, сводится к задаче вычисления значения некоторой функции 

УДК 519.65

М.А. МАЙДАКОВ

M.A. Maidakov

Дискретная интерполяция методом естественных соседей на основе подхода с разбросом значений

Discrete Natural Neighbor Interpolation method based on scatter approach

Аннотация: В работе рассмотрен модифицированный метод построения дискретной интерполяции естественных соседей на основе подхода с «разбросом» значения, однако не требующий предварительного построения деревьев для обхода и поиска ближайшего соседа. Основная идея предлагаемого метода заключается в объединении задач поиска значения ближайшего соседа и разброса значения. Данный метод хорошо подходит для распараллеливания вычислений на системах с общей памятью, в том числе и на графических процессорах с поддержкой технологий CUDA и OpenCL.

Ключевые слова: дискретная интерполяция естественных соседей, интерполяция Сибсона, диаграмма Вороного, дискретное представление, параллельные вычисления на графических процессорах.

Summary: This paper describes a modified method for constructing discrete natural neighbors interpolation based on the "scatter" approach, but does not require pre-construction trees for searching the nearest neighbor. The main idea of the proposed method is to combine the tasks of finding the values of the nearest neighbor and scatter values. This method is well suited for parallel computing systems with shared memory, including graphics processors with support for CUDA or OpenCL.
Keywords: discrete natural neighbors interpolation, Sibson interpolation, Voronoi diagram, discrete representation, parallel computing on GPU.

Задача интерполяции, в общем случае, сводится к задаче вычисления значения некоторой функции ₀ = (x₀) в заданной точке x₀, по ее значениям, заданным на фиксированной системе узлов x_i  E^d. Алгоритмы вычисления значений неизвестной функции  = (x), тем или иным способом сводятся к окончательной формуле вида:

, где :

_i = 1 (1)

В выражении (1) a_i – коэффициенты интерполяции, зависящие от расположения системы точек и не зависящие от самой функции ; n – количество точек, по которым производится интерполяция функции .

Различные алгоритмы интерполяции отличаются способом выбора коэффициентов интерполяции и методами выбора соседей для узловой точки[1]. Интерполяция естественных соседей существенного лучше традиционных методов дистанционного взвешивания, так как её плотность меняется в зависимости от суммарного вклада влияющих областей, а также обеспечивается единственность и непрерывность результатов интерполяции всюду за исключением мест опробования.

Интерполяция по методу естественных соседей была введена Сибсоном в 1981 году как средство для приближения и сглаживания рассеянных данных. В общем случае интерполянт Сибсона функции  оценивается в узле p следующим образом [2]:

, где

(2)

Где u_i – является пересечением площадей (объемов) ячейки Воронного для точки p_i и временной ячейки Вороного для точки p.

При этом диаграмма Вороного представляет собой разбиение некоторой области на ячейки или регионы на основе заданного конечного множества рассеянных точек данных, названных сайтами. Диаграмма Вороного U(N) множества N в области  является разбиением области на регионы U(p_i)   таким образом, что любая точка в U(p_i) находится ближе к пробе pi, чем к любой другой пробе p_j  N(j  i). Регион U(p_i) ассоциированный с пробой p_i называется ячейкой Вороного и определяется как:

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

Традиционные подходы построения дискретной интерполяции Сибсона основаны на вычислении весовых вкладов интерполянта путём определения площадей (в двумерном случае) или объемов (в трехмерном случае) u_i, связанных с пробами p_i после вставки в диаграмму Вороного временного узла p[3]. Поскольку расчет интерполянта Сибсона в узле p эквивалентен определению средних значений соседних проб, то дискретная интерполяция Сибсона может быть вычислена путём усреднения дискретных элементов. Накапливая все значения данных опробования для исходной диаграммы Вороного внутри временно вставляемого региона Вороного U(p) и разделив накопленное значение на число просуммированных элементов, мы определим среднее значение для региона. Это значение и есть дискретный эквивалент значения интерполяции Сибсона. Поскольку дискретные элементы накапливаются в одной позиции, то данная схема интерполяции называется “подходом собирания”.

Более эффективный подход к построению дискретной интерполяции по методу естественных соседей был предложен в работе [3]. Ключевая идея построения состоит в рассмотрении проблемы с “обратной” стороны, когда значения “разбрасываются”, а не собираются. Тогда вместо итеративного прохода по всем позициям p и накопления всех значений позиций растра i  V(p) в c(p), можно выполнить итерации по растровым позициям i и определять, какие (p) находятся под влиянием значения растра i. Смешав d-мерные сферы для каждого растра i в результате мы получим дискретный интерполянт Сибсона. Нахождение ближайшего соседа происходит без прямого построения диаграммы Воронго, а при помощи построения kd-дерева.

Данный подход намного быстрее и проще в реализации по сравнению с традиционными подходами, поскольку не требует вычисления диаграммы Вороного для каждого растра (вокселя) регулярной решетки. Однако построение бинарных деревьев для поиска ближайшего соседа все еще остается слабым местом данного подхода. Для построение kd-дерева для n точек на предварительном этапе требуется 0(n*logn) времени. Запрос ближайшего соседа для каждого растра i потребует еще O(logn) времени. При довольно больших наборах исходных данных, время поиска ближайшего соседа каждого растра (вокселя) будет слишком велико.

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

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

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

Разброс указателя на значение для каждого растра осуществляется по поверхности n-мерной сферы с последовательным увеличением радиуса разброса, и на каждой итерации необходимо вычисление позиции разбрасываемого значения на n-мерной сфере. Тогда для каждого растра (вокселя) регулярной решетки будет вычислено πR² позиций для двумерного случая и 4πR² для трехмерного, где R – конечный радиус для каждого растра (вокселя).

Поскольку для каждого узла регулярной решетки разброс значений осуществляется одинаковы образом, то вычисление позиций для разброса можно вычислить на начальном этапе и представить в виде последовательности смещений от узла регулярной решетки. При этом можно сократить количество вычислений для поиска позиций разброса в N*M раз для двумерного случая и N*M*L для трехмерного, где N,M,L –размеры регулярной решетки.

Тогда алгоритм построения дискретной интерполяции по методу естественных соседей для CPU будет выглядеть следующим образом:

Инициализация входных данных.
Заполнение массива RP с последовательностью разброса.
Для каждого узла p регулярной решетки:
1. Пока не найден ближайший сосед:
  - Разбрасывать указатель на p для каждого узла в соответствии с последовательностью RP.
2. Записать значение найденного ближайшего соседа в ячейку памяти p.
Для каждого узла p регулярной решетки:
1. Получить усредненное значение для каждой ячейки равное делению суммы накопленных значений на количество накоплений.
Вывод результатов интерполяции.

Для вычислительных систем на основе GPU из-за невозможности динамического выделения памяти в процессе выполнения предлагается заменить: разброс указателя, разбросом значения на обратном проходе.

Тогда для вычислительных систем на основе GPU необходимо заменить пункт 3 в описанном выше алгоритме следующим пунктом:

3. Для каждого узла p регулярной решетки:

a) Пока не найден ближайший сосед:

Установить счетчик с = 0.
Проверить каждый узел в соответствии с последовательностью RP на наличие значения.
Увеличить счетчик с на единицу.

b) Пока счетчик с >= 0:

Добавить найденное значение ближайшего соседа к значению в ячейках в соответствии с последовательностью разброса RP[c].
Уменьшить с на единицу.

Предложенный метод, в отличие от стандартных методов построения дискретной интерполяции естественных соседей, достаточно прост в реализации и не использует сложных структур данных для хранения промежуточных вычислений, поскольку нет необходимости в построении деревьев для поиска естественных соседей, а также в вычислении площадей (объемов) пересечения ячеек Вороного. При этом данный метод достаточно хорошо подходит для распараллеливания вычислений на системах с общей памятью, так как этапы 3 и 4 в приведенном выше алгоритме могут выполняться для каждого узла регулярной решетки независимо, и скорость вычисления будет в большей степени зависеть от количества одновременно выполняющихся процессов.

СПИСОК ЛИТЕРАТУРЫ

Васильев, П.В. Использование графического ускорения интерполяции Сибсона для моделирования геотекстур [Текст] / П.В. Васильев, Майдаков М.А. // I-я международная научно-техническая конференция «Компьютерные науки и технологии». Издательство Белгородского государственного университета. -2009. –C.67-69.
Sibson, R. A brief description of natural neighbor interpolation in Multivariate Data / R. Sibson // Jhon Wiley & Sons, 1981. - P. 21-36.
Sung W. Park. Discrete Sibson Interpolation./ Sung W. Park, Lars Linsen, Oliver Kreylos, John D. Owens, Bernd Hamann.// IEEE Transactions On Visualization And Computer Graphics, Vol. 12, No. 2, 2006.
Liang L, Hale D. A stable and fast implementation of natural neighbor interpolation / Liang L, Hale D. // Center for Wave Phenomena, Colorado School of Mines,Golden, CO 80401, USA, -2010. –P. 207-219.

Майдаков Михаил Александрович

Белгородский государственный университет, г. Белгород

Аспирант кафедры математического и программного обеспечения

информационных систем

Тел.: 8-919-2225525

E-mail: m-email@inbox.com