I've tried to replicate the main bottleneck in one of my programs.
I want to get the linearly (or rather bilinearly) interpolated values of several non-integer pixel values simultaneously. It is not the case that each pixel coordinate is perturbed in the same way. Below is a complete/minimal script along with comments that demonstrates the problem. How can I speed up the calculation of result
?
import numpy as np
import time
im = np.random.rand(640,480,3) # my "image"
xx, yy = np.meshgrid(np.arange(im.shape[1]), np.arange(im.shape[0]))
print "Check these are the right indices:",np.sum(im - im[yy,xx,:])
# perturb the indices slightly
# I want to calculate the interpolated
# values of "im" at these locations
xx = xx + np.random.normal(size=im.shape[:2])
yy = yy + np.random.normal(size=im.shape[:2])
# integer value/pixel locations
x_0 = np.int_(np.modf(xx)[1])
y_0 = np.int_(np.modf(yy)[1])
x_1, y_1 = x_0 + 1, y_0 + 1
# the real-valued offsets/coefficients pixels
a = np.modf(xx)[0][:,:,np.newaxis]
b = np.modf(yy)[0][:,:,np.newaxis]
# make sure we don't go out of bounds at edge pixels
np.clip(x_0,0,im.shape[1]-1,out=x_0)
np.clip(x_1,0,im.shape[1]-1,out=x_1)
np.clip(y_0,0,im.shape[0]-1,out=y_0)
np.clip(y_1,0,im.shape[0]-1,out=y_1)
# now perform linear interpolation: THIS IS THE BOTTLENECK!
tic = time.time()
result = ((1-a) * (1-b) * im[y_0, x_0, :] +
a * (1-b) * im[y_1, x_0, :] +
(1-a) * b * im[y_0, x_1, :] +
a * b * im[y_1, x_1, :] )
toc = time.time()
print "interpolation time:",toc-tic
The simplest solution is to use something which can be vectorized. eg. ynew = function(xnew);simps(ynew,xnew) This is much faster, but depending on the inputs less accurate. Another possibility which is also a lot faster and gives the same results is to implement a low level callable.
interp() function returns the one-dimensional piecewise linear interpolant to a function with given discrete data points (xp, fp), evaluated at x. Parameters : x : [array_like] The x-coordinates at which to evaluate the interpolated values.
Linear Interpolation is the technique of determining the values of the functions of any intermediate points when the values of two adjacent points are known. Linear interpolation is basically the estimation of an unknown value that falls within two known values.
Linear Interpolation means estimating a missing value by connecting dots in the straight line in increasing order. It estimates the unknown value in the same increasing order as the previous values. The default method used by Interpolation is Linear so while applying it one does not need to specify it.
Thanks to @JoeKington for the suggestion. Here's the best I can come up with using scipy.ndimage.map_coordinates
# rest as before
from scipy import ndimage
tic = time.time()
new_result = np.zeros(im.shape)
coords = np.array([yy,xx,np.zeros(im.shape[:2])])
for d in range(im.shape[2]):
new_result[:,:,d] = ndimage.map_coordinates(im,coords,order=1)
coords[2] += 1
toc = time.time()
print "interpolation time:",toc-tic
Update: Added the tweaks suggested in the comments and tried one or two other things. This is the fastest version:
tic = time.time()
new_result = np.zeros(im.shape)
coords = np.array([yy,xx])
for d in range(im.shape[2]):
ndimage.map_coordinates(im[:,:,d],
coords,order=1,
prefilter=False,
output=new_result[:,:,d] )
toc = time.time()
print "interpolation time:",toc-tic
Example running time:
original version: 0.463063955307
better version: 0.204537153244
best version: 0.121845006943
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With