Cython.pyx

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
from __future__ import division

import numpy as np
cimport numpy as np

from libc.stdlib cimport malloc, free

from math import sqrt


cdef float cerror_pixels(long[:,:] pixels, float[:,:] distances, int dx, int dy):
    cdef int w = distances.shape[0]
    cdef int h = distances.shape[1]
    cdef int n = pixels.shape[0]

    cdef float dsum = 0

    cdef int k = 0
    cdef int x = 0
    cdef int y = 0
    cdef int n_inbound = 0

    cdef float* pixel_distances = <float*> malloc(sizeof(float)*n)
    cdef float float_tmp = 0

    for k in range(n):
        x = pixels[k,0] + dx
        y = pixels[k,1] + dy
        if (x >= 0) and (y >= 0) and (x < w) and (y < h):
            float_tmp = distances[x,y]
            float_tmp *= float_tmp

            pixel_distances[n_inbound] = float_tmp
            dsum += float_tmp
            n_inbound += 1
    
    cdef float result = 0

    cdef float dmean = 0
    cdef float dstddev = 0
    cdef float dstderr = 0

    if n_inbound > 0:
        dmean = (<float>dsum) / (<float>n_inbound)
        dstddev = 0
        for k in range(n_inbound):
            float_tmp = (pixel_distances[k] - dmean)
            dstddev += float_tmp*float_tmp
        dstddev = sqrt(dstddev/n_inbound)
        dstderr = dstddev / sqrt(n_inbound)
        result = dmean + 2*dstderr
    else:
        result = 1e5 # big value but not too big or we lose precision
    
    free(pixel_distances)
    return result

def cerror_pixels_wrapper(pixels,distances,dx=0,dy=0):
    return cerror_pixels(pixels,distances,dx,dy)

Python.py

def error_pixels(pixels,distances):
    w,h = distances.shape[:2]
    in_bound = np.logical_and(
                    np.logical_and(pixels[:,0] >= 0,pixels[:,1] >= 0),
                    np.logical_and(pixels[:,0] < w,pixels[:,1] < h))
    ds = distances[pixels[in_bound,0],pixels[in_bound,1]]
    ds2 = ds*ds
    dsum = np.sum(ds2)
    
    n = ds.shape[0]
    if n > 0:
        dstd = np.std(ds2)
        dmean = dsum/n
        dstderr = dstd/sqrt(n)
        return dmean+2*dstderr
    else:
        return 1e5 # big value but not too big
        #return dsum, n

Speed up von gut 20x:

In [62]:
%timeit cerror_pixels_wrapper(_pixels,dist["distances"])
100000 loops, best of 3: 4.34 µs per loop

In [63]:
%timeit error_pixels(_pixels,dist["distances"])
10000 loops, best of 3: 99.7 µs per loop

In [64]:
99.7/4.34
Out[64]:
22.972350230414747