Given a point (0.8, 0.6) and an intensity (3), one can bi-linearly “reverse interpolate” a single point onto a 2×2 grid with integer indices (0,0) -> (1,1).
The points and grid are oriented with the first entry increasing downwards, and the second entry increasing towards the right.
On such a grid, the weight of just the coordinate above becomes:
0.08 | 0.12 ---------------- 0.32 | 0.48
and we can multiply by the intensity associated with the coordinate to give a 2×2 grid bilinearly weighted intensity:
0.24 | 0.36 ---------------- 0.96 | 1.44
And can be plotted like this:
For several points, one can weigh these onto the same array (full code below):
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)
For my work, I require both the weights and the intensity*weights as output arrays. I need to perform the above on a very large (10s of millions) number of coordinates, where the coordinates have values from 0.0 to 4095.0. I have written a numpy routine below, and while it is reasonably fast (1.25 sec) for 100_000 points, I would prefer it was faster as I need to call it several times on the 10_000_000 data points I have.
I considered vectorising the numpy code instead of a for loop, but then I’m generating a 4096×4096 mostly empty array for each point that I would then sum. That would require 1000 TB of ram.
I have also tried a naive implementation in cupy, but since I employ a for loop, it becomes far too slow.
In my code, I generate a 2×2 weighted array for each point, multiply the array by the intensity, and add those to the main arrays by slicing. Is there a better way?
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)
For testing the code, I also offer the following plotting code – only use with a small (~10) number of points, or the scatter will cover the whole image.
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')
Advertisement
Answer
Thanks @armamut for a good answer! It inspired me to look a bit, and I found np.bincount, which is also implemented in cupy. It turns out that the bincount implementation is faster, and the cupy implementation is really fast! The latter could probably be improved a bit further, since I had to juggle a few tuples to get it to work.
# Timings points = np.random.random((10_000_000, 2)) * (256, 256) intens = np.random.random((10_000_000)) pcupy = cp.asarray(points) icupy = cp.asarray(intens) %time bilinear_bincount_cupy(pcupy, icupy) %time bilinear_bincount_numpy(points, intens) %time bilinear_2(points, intens) Wall time: 456 ms Wall time: 2.57 s Wall time: 5.37 s
The numpy implementation:
def bilinear_bincount_numpy(points, intensities):
"""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)
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)
shifts = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
indices = floored_indices[:, None] + shifts
indices = (indices * (shape[1], 1)).sum(-1)
weights = np.array([w1, w2, w3, w4]).T
weight_bins = np.bincount(indices.flatten(), weights=weights.flatten())
intens_bins = np.bincount(indices.flatten(), weights=(intensities[:, None]*weights).flatten())
all_weight_bins = np.zeros(np.prod(shape))
all_intens_bins = np.zeros(np.prod(shape))
all_weight_bins[:len(weight_bins)] = weight_bins
all_intens_bins[:len(weight_bins)] = intens_bins
weight_image = all_weight_bins.reshape(shape)
intens_image = all_intens_bins.reshape(shape)
return intens_image, weight_image
And the cupy implementation:
def bilinear_bincount_cupy(points, intensities):
"""Bilinear weighting of points onto a grid.
Extent of grid given by min and max of points in each dimension
points should be a cupy array of shape (N, 2)
intensity should be a cupy array of shape (N,)
"""
floor = cp.floor(points)
ceil = floor + 1
floored_indices = cp.array(floor, dtype=int)
low0, low1 = floored_indices.min(0)
high0, high1 = floored_indices.max(0)
floored_indices = floored_indices - cp.array([low0, low1])
shape = cp.array([high0 - low0 + 2, high1-low1 + 2])
upper_diff = ceil - points
lower_diff = points - floor
w1 = upper_diff[:, 0] * upper_diff[:, 1]
w2 = upper_diff[:, 0] * lower_diff[:, 1]
w3 = lower_diff[:, 0] * upper_diff[:, 1]
w4 = lower_diff[:, 0] * lower_diff[:, 1]
shifts = cp.array([[0, 0], [0, 1], [1, 0], [1, 1]])
indices = floored_indices[:, None] + shifts
indices = (indices * cp.array([shape[1].item(), 1])).sum(-1)
weights = cp.array([w1, w2, w3, w4]).T
# These bins only fill up to the highest index - not to shape[0]*shape[1]
weight_bins = cp.bincount(indices.flatten(), weights=weights.flatten())
intens_bins = cp.bincount(indices.flatten(), weights=(intensities[:, None]*weights).flatten())
# So we create a zeros array that is big enough
all_weight_bins = cp.zeros(cp.prod(shape).item())
all_intens_bins = cp.zeros_like(all_weight_bins)
# And fill it here
all_weight_bins[:len(weight_bins)] = weight_bins
all_intens_bins[:len(weight_bins)] = intens_bins
weight_image = all_weight_bins.reshape(shape.get())
intens_image = all_intens_bins.reshape(shape.get())
return intens_image, weight_image