benchmarks.py

import numpy as np
from lib import get_scale_local_maximas_orignal, get_scale_local_maximas_vectorized
from compiled import get_scale_local_maximas_cython

lapl_dummy = np.random.rand(100,100,100)
coords = np.random.random_integers(0,99, size=(1000,3))

%timeit get_scale_local_maximas_orignal(coords, lapl_dummy)
%timeit get_scale_local_maximas_cython(coords, lapl_dummy)

%timeit get_scale_local_maximas_vectorized(coords, lapl_dummy)

# Output:
#1000 loops, best of 3: 1.6 ms per loop
#1000 loops, best of 3: 328 µs per loop
#10000 loops, best of 3: 103 µs per loop

compiled.pyx

# cython: cdivision=True
# cython: boundscheck=False
# cython: nonecheck=False
# cython: wraparound=False
import numpy as np
cimport numpy as cnp


def get_scale_local_maximas_cython(cnp.ndarray[cnp.int_t, ndim=2] cube_coordinates, cnp.ndarray[cnp.double_t, ndim=3] laplacian_cube):
    """
    Check provided cube coordinate for scale space local maximas.
    Returns only the points that satisfy the criteria.

    A point is considered to be a local maxima if its value is greater
    than the value of the point on the next scale level and the point
    on the previous scale level. If the tested point is located on the
    first scale level or on the last one, then only one inequality should
    hold in order for this point to be local scale maxima.

    Parameters
    ----------
    cube_coordinates : (n, 3) ndarray
          A 2d array with each row representing 3 values, ``(y,x,scale_level)``
          where ``(y,x)`` are coordinates of the blob and ``scale_level`` is the
          position of a point in scale space.
    laplacian_cube : ndarray of floats
        Laplacian of Gaussian scale space. 

    Returns
    -------
    output : (n, 3) ndarray
        cube_coordinates that satisfy the local maximum criteria in
        scale space.

    Examples
    --------
    >>> one = np.array([[1, 2, 3], [4, 5, 6]])
    >>> two = np.array([[7, 8, 9], [10, 11, 12]])
    >>> three = np.array([[0, 0, 0], [0, 0, 0]])
    >>> check_coords = np.array([[1, 0, 1], [1, 0, 0], [1, 0, 2]])
    >>> lapl_dummy = np.dstack([one, two, three])
    >>> get_scale_local_maximas(check_coords, lapl_dummy)
    array([[1, 0, 1]])
    """

    cdef Py_ssize_t y_coord, x_coord, point_layer, point_index
    cdef cnp.double_t point_response, lower_point_response, upper_point_response
    cdef Py_ssize_t amount_of_layers = laplacian_cube.shape[2]
    cdef Py_ssize_t amount_of_points = cube_coordinates.shape[0]

    # Preallocate index. Fill it with False.
    accepted_points_index = np.ones(amount_of_points, dtype=bool)

    for point_index in range(amount_of_points):

        interest_point_coords = cube_coordinates[point_index]
        # Row coordinate
        y_coord = interest_point_coords[0]
        # Column coordinate
        x_coord = interest_point_coords[1]
        # Layer number starting from the smallest sigma
        point_layer = interest_point_coords[2]
        point_response = laplacian_cube[y_coord, x_coord, point_layer]

        # Check the point under the current one
        if point_layer != 0:
            lower_point_response = laplacian_cube[y_coord, x_coord, point_layer-1]
            if lower_point_response >= point_response:
                accepted_points_index[point_index] = False
                continue

        # Check the point above the current one
        if point_layer != (amount_of_layers-1):
            upper_point_response = laplacian_cube[y_coord, x_coord, point_layer+1]
            if upper_point_response >= point_response:
                accepted_points_index[point_index] = False
                continue

    # Return only accepted points
    return cube_coordinates[accepted_points_index]

lib.py

def get_scale_local_maximas_orignal(cube_coordinates, laplacian_cube):
    """
    Check provided cube coordinate for scale space local maximas.
    Returns only the points that satisfy the criteria.

    A point is considered to be a local maxima if its value is greater
    than the value of the point on the next scale level and the point
    on the previous scale level. If the tested point is located on the
    first scale level or on the last one, then only one inequality should
    hold in order for this point to be local scale maxima.

    Parameters
    ----------
    cube_coordinates : (n, 3) ndarray
          A 2d array with each row representing 3 values, ``(y,x,scale_level)``
          where ``(y,x)`` are coordinates of the blob and ``scale_level`` is the
          position of a point in scale space.
    laplacian_cube : ndarray of floats
        Laplacian of Gaussian scale space. 

    Returns
    -------
    output : (n, 3) ndarray
        cube_coordinates that satisfy the local maximum criteria in
        scale space.

    Examples
    --------
    >>> one = np.array([[1, 2, 3], [4, 5, 6]])
    >>> two = np.array([[7, 8, 9], [10, 11, 12]])
    >>> three = np.array([[0, 0, 0], [0, 0, 0]])
    >>> check_coords = np.array([[1, 0, 1], [1, 0, 0], [1, 0, 2]])
    >>> lapl_dummy = np.dstack([one, two, three])
    >>> get_scale_local_maximas(check_coords, lapl_dummy)
    array([[1, 0, 1]])
    """

    amount_of_layers = laplacian_cube.shape[2]
    amount_of_points = cube_coordinates.shape[0]

    # Preallocate index. Fill it with False.
    accepted_points_index = np.ones(amount_of_points, dtype=bool)

    for point_index, interest_point_coords in enumerate(cube_coordinates):
        # Row coordinate
        y_coord = interest_point_coords[0]
        # Column coordinate
        x_coord = interest_point_coords[1]
        # Layer number starting from the smallest sigma
        point_layer = interest_point_coords[2]
        point_response = laplacian_cube[y_coord, x_coord, point_layer]

        # Check the point under the current one
        if point_layer != 0:
            lower_point_response = laplacian_cube[y_coord, x_coord, point_layer-1]
            if lower_point_response >= point_response:
                accepted_points_index[point_index] = False
                continue

        # Check the point above the current one
        if point_layer != (amount_of_layers-1):
            upper_point_response = laplacian_cube[y_coord, x_coord, point_layer+1]
            if upper_point_response >= point_response:
                accepted_points_index[point_index] = False
                continue

    # Return only accepted points
    return cube_coordinates[accepted_points_index]

def get_scale_local_maximas_vectorized(cube_coordinates, laplacian_cube):
    """
    Check provided cube coordinate for scale space local maximas.
    Returns only the points that satisfy the criteria.

    A point is considered to be a local maxima if its value is greater
    than the value of the point on the next scale level and the point
    on the previous scale level. If the tested point is located on the
    first scale level or on the last one, then only one inequality should
    hold in order for this point to be local scale maxima.

    Parameters
    ----------
    cube_coordinates : (n, 3) ndarray
          A 2d array with each row representing 3 values, ``(y,x,scale_level)``
          where ``(y,x)`` are coordinates of the blob and ``scale_level`` is the
          position of a point in scale space.
    laplacian_cube : ndarray of floats
        Laplacian of Gaussian scale space. 

    Returns
    -------
    output : (n, 3) ndarray
        cube_coordinates that satisfy the local maximum criteria in
        scale space.

    Examples
    --------
    >>> one = np.array([[1, 2, 3], [4, 5, 6]])
    >>> two = np.array([[7, 8, 9], [10, 11, 12]])
    >>> three = np.array([[0, 0, 0], [0, 0, 0]])
    >>> check_coords = np.array([[1, 0, 1], [1, 0, 0], [1, 0, 2]])
    >>> lapl_dummy = np.dstack([one, two, three])
    >>> get_scale_local_maximas(check_coords, lapl_dummy)
    array([[1, 0, 1]])
    """
    x, y, z = [ cube_coordinates[:, ind] for ind in range(3) ]

    point_responses = laplacian_cube[x, y, z]
    lowers = point_responses.copy()
    uppers = point_responses.copy()
    not_layer_0 = z > 0
    lower_responses = laplacian_cube[x[not_layer_0], y[not_layer_0], z[not_layer_0]-1]
    lowers[not_layer_0] = lower_responses  

    not_max_layer = z < (laplacian_cube.shape[2] - 1)
    upper_responses = laplacian_cube[x[not_max_layer], y[not_max_layer], z[not_max_layer]+1]
    uppers[not_max_layer] = upper_responses

    lo_check = np.ones(z.shape, dtype=np.bool)
    lo_check[not_layer_0] = (point_responses > lowers)[not_layer_0]
    hi_check = np.ones(z.shape, dtype=np.bool)
    hi_check[not_max_layer] = (point_responses > uppers)[not_max_layer]

    return cube_coordinates[lo_check & hi_check]