CODEWITHSUNDEEP
pythonc++numpyfilter
Alexander McFarlane
Alexander McFarlane

Reputation: 11293

A faster discrete Laplacian than scipy.ndimage.filters.laplace for small arrays

My spends the vast bulk of its computational time in scipy.ndimage.filters.laplace()

The main advantage of scipy and numpy is vectorised calculation in C/C++ wrapped in python. scipy.ndimage.filters.laplace() is derived from _nd_image.correlate1d which is part of the optimised library nd_image.h

Is there any faster method of doing this across an array of size 10-100?

Definition Laplace Filter - ignoring division

Upvotes: 2

Views: 1651

Answers (1)

Alexander McFarlane
Alexander McFarlane

Reputation: 11293

The problem was rooted in scipy's excellent error handling and debugging. However, in the instance the user knows what they're doing it just provides excess overhead.

This code below strips all the python clutter in the back end of scipy and directly accesses the C++ function to get a ~6x speed up!

laplace == Mine ? True
testing timings...
array size 10
100000 loops, best of 3: 12.7 µs per loop
100000 loops, best of 3: 2.3 µs per loop
array size 100
100000 loops, best of 3: 12.7 µs per loop
100000 loops, best of 3: 2.5 µs per loop
array size 100000
1000 loops, best of 3: 413 µs per loop
1000 loops, best of 3: 404 µs per loop

Code

from scipy import ndimage
from scipy.ndimage import _nd_image
import numpy as np

laplace_filter = np.asarray([1, -2, 1], dtype=np.float64)

def fastLaplaceNd(arr):
    output = np.zeros(arr.shape, 'float64')
    if arr.ndim > 0:
        _nd_image.correlate1d(arr, laplace_filter, 0, output, 1, 0.0, 0)
        if arr.ndim == 1: return output
        for ax in xrange(1, arr.ndim):
            output += _nd_image.correlate1d(arr, laplace_filter, ax, output, 1, 0.0, 0)
    return output

if __name__ == '__main__':
    arr = np.random.random(10)
    test = (ndimage.filters.laplace(arr, mode='wrap') == fastLaplace(arr)).all()
    assert test
    print "laplace == Mine ?", test
    print 'testing timings...'
    print "array size 10"
    %timeit ndimage.filters.laplace(arr, mode='wrap')
    %timeit fastLaplace(arr)
    print 'array size 100'
    arr = np.random.random(100)
    %timeit ndimage.filters.laplace(arr, mode='wrap')
    %timeit fastLaplace(arr)
    print "array size 100000"
    arr = np.random.random(100000)
    %timeit ndimage.filters.laplace(arr, mode='wrap')
    %timeit fastLaplace(arr)

Upvotes: 2

Related Questions