Имея точку (0.8, 0.6) и интенсивность (3), можно выполнить билинейную обратную интерполяцию одной точки на сетку 2x2 с целочисленными индексами (0,0) -> (1,1).
Точки и сетка ориентированы так, что первая запись увеличивается вниз, а вторая запись увеличивается вправо.
В такой сетке вес только указанной выше координаты становится:
0.08 | 0.12
----------------
0.32 | 0.48
и мы можем умножить на интенсивность, связанную с координатой, чтобы получить билинейно взвешенную интенсивность сетки 2x2:
0.24 | 0.36
----------------
0.96 | 1.44
А можно изобразить так:
Для нескольких точек можно взвесить их в одном и том же массиве (полный код ниже):
points = np.array([(0.8, 0.6), (2.2, 2.6),(5, 1), (3, 4.2), (8.5, 8.2)])
intens = np.array([3, 3, 1, 1, 2])
image, weight = bilinear(points, intens)
Для моей работы мне нужны как weights, так и intensity*weights в качестве выходных массивов. Мне нужно выполнить вышеуказанное для очень большого (десятки миллионов) числа координат, где координаты имеют значения от 0,0 до 4095,0. Я написал подпрограмму numpy ниже, и хотя она достаточно быстра (1,25 секунды) для 100_000 точек, я бы предпочел, чтобы она была быстрее, так как мне нужно вызывать ее несколько раз для 10_000_000 точек данных, которые у меня есть.
Я подумал о векторизации кода numpy вместо цикла for, но затем я генерирую в основном пустой массив 4096x4096 для каждой точки, которую затем суммирую. Для этого потребуется 1000 ТБ оперативной памяти.
Я также пробовал наивную реализацию в cupy, но поскольку я использую цикл for, он становится слишком медленным.
В моем коде я генерирую взвешенный массив 2x2 для каждой точки, умножаю массив на интенсивность и добавляю их к основным массивам путем нарезки. Есть ли способ лучше?
import numpy as np
def bilinear(points, intensity):
"""Bilinear weighting of points onto a grid.
Extent of grid given by min and max of points in each dimension
points should have shape (N, 2)
intensity should have shape (N,)
"""
floor = np.floor(points)
ceil = floor + 1
floored_indices = np.array(floor, dtype=int)
low0, low1 = floored_indices.min(0)
high0, high1 = floored_indices.max(0)
floored_indices = floored_indices - (low0, low1)
shape = (high0 - low0 + 2, high1-low1 + 2)
weights_arr = np.zeros(shape, dtype=float)
int_arr = np.zeros(shape, dtype=float)
upper_diff = ceil - points
lower_diff = points - floor
w1 = np.prod((upper_diff), axis=1)
w2 = upper_diff[:,0]*lower_diff[:,1]
w3 = lower_diff[:,0]*upper_diff[:,1]
w4 = np.prod((lower_diff), axis=1)
for i, index in enumerate(floored_indices):
s = np.s_[index[0]:index[0]+2, index[1]:index[1]+2]
weights = np.array([[w1[i], w2[i]], [w3[i], w4[i]]])
weights_arr[s] += weights
int_arr[s] += intensity[i]*weights
return int_arr, weights_arr
rng = np.random.default_rng()
N_points = 10_000 # use 10_000 so it is quick
image_shape = (256, 256) # Use 256 so it isn't so big
points = rng.random((N_points, 2)) * image_shape
intensity = rng.random(N_points)
image, weight = bilinear(points, intensity)
Для тестирования кода я также предлагаю следующий код построения графика — используйте только с небольшим (~ 10) количеством точек, иначе разброс покроет все изображение.
import matplotlib.pyplot as plt
floor = np.floor(points) - 0.5
lower, left = floor.min(0)
upper, right = (floor).max(0) + 2
extent = (left, right, upper, lower)
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(6,3))
ax1.scatter(*points[:,::-1].T, c='red')
im1 = ax1.imshow(weight, clim=(image.min(), image.max()), extent=extent)
ax1.set(title='Weight', xlim=(left - 1, right + 1), ylim = (upper + 1, lower - 1))
colorbar(im1)
ax2.scatter(*points[:,::-1].T , c='red')
im2 = ax2.imshow(image, extent=extent)
ax2.set(title='Weight x Intensity', xlim=(left - 1, right + 1), ylim = (upper + 1, lower - 1))
colorbar(im2)
plt.tight_layout()
plt.show()
# If labeling the first point
# ax1.text(*points[0].T, f"({points[0,0]}, {points[0,1]})", va='bottom', ha='center', color='red')
# ax2.text(*points[0].T, f"({points[0,0]}, {points[0,1]}, {intens[0]})", va='bottom', ha='center', color='red')
