Как использовать scipy.interpolate.LinearNDInterpolator с собственной триангуляцией

Question

У меня есть собственный алгоритм триангуляции, который создает триангуляцию, основанную как на условии Делоне, так и на градиенте, так что треугольники совпадают с градиентом.

Это пример вывода:

Приведенное выше описание не имеет отношения к вопросу, но необходимо для контекста.

Теперь я хочу использовать мою триангуляцию с scipy.interpolate.LinearNDInterpolator сделать интерполяцию.

С Делони Сципи я бы сделал следующее

import numpy as np
import scipy.interpolate
import scipy.spatial
points = np.random.rand(100, 2)
values = np.random.rand(100)
delaunay = scipy.spatial.Delaunay(points)
ip = scipy.interpolate.LinearNDInterpolator(delaunay, values)

это delaunay объект имеет delaunay.points а также delaunay.simplices которые образуют триангуляцию. У меня точно такая же информация с моей собственной триангуляцией, но scipy.interpolate.LinearNDInterpolator требует scipy.spatial.Delaunay объект.

Я думаю, что мне нужно было бы подкласс scipy.spatial.Delaunay и реализовать соответствующие методы. Однако я не знаю, какие мне нужны, чтобы туда попасть.

Answer 1

Я хотел сделать то же самое с триангуляцией Делоне, предложеннойtriangleпакет . Треугольный код Делоне примерно в восемь раз быстрее кода SciPy на больших (~100_000) точках. (Я призываю других разработчиков попытаться превзойти это :))

К сожалению,Scipy LinearNDInterpolatorФункция в значительной степени зависит от определенных атрибутов, присутствующих в объекте триангуляции Делоне SciPy. Они создаются_get_delaunay_info()Код CPython, который трудно дизассемблировать. Даже зная, какие атрибуты необходимы (похоже, их много, включая такие вещи, как paraboloid_scaleа также paraboloid_shift), я не знаю, как извлечь это из другой библиотеки триангуляции.

Вместо этого я попробовал подход @Patol75 (комментарий к вопросу выше), но используяLinearTriInterpolatorвместо Кубика. Код работает правильно, но медленнее, чем в SciPy. Интерполяция 400_000 точек из облака 400_000 точек занимает примерно в 3 раза больше времени с использованием кода matplotlib, чем scipy. Код Matplotlib tri написан на C++, поэтому преобразование кода в CuPy не так просто. Если бы мы могли смешать два подхода, мы могли бы сократить общее время с 3,65 с / 10,2 с до 1,1 секунды!

      import numpy as np
np.random.seed(1)
N = 400_000
shape = (100, 100)
points = np.random.random((N, 2)) * shape # spread over 100, 100 to avoid float point errors
vals = np.random.random((N,))
interp_points1 = np.random.random((N,2)) * shape
interp_points2 = np.random.random((N,2)) * shape

triangle_input = dict(vertices=points)

### Matplotlib Tri
import triangle as tr
from matplotlib.tri import Triangulation, LinearTriInterpolator

triangle_output = tr.triangulate(triangle_input) # 280 ms
tri = tr.triangulate(triangle_input)['triangles'] # 280 ms
tri = Triangulation(*points.T, tri) # 5 ms
func = LinearTriInterpolator(tri, vals) # 9490 ms
func(*interp_points.T).data # 116 ms
# returns [0.54467719, 0.35885304, ...]
# total time 10.2 sec

### Scipy interpolate
tri = Delaunay(points) # 2720 ms
func = LinearNDInterpolator(tri, vals) # 1 ms
func(interp_points) # 925 ms

# returns [0.54467719, 0.35885304, ...]
# total time 3.65 sec