I coded a bilateral filter class in cython based on numpy's neighborhood
iterator (thanks to T. J's code example). I tried to parallel the code by
replacing the standard loop (commented line 150) by a prange loop (line 151).
The result are series of compilation errors mainly due the the use of
iterators. Is there an *easy* way to work with nmpy iterators while the GIL is
released?
Platfoem: numpy 1.6.1 on python 2.7.2 and cython 0.15.1
system: gcc on linux
Nadav
# A Cython + Neighbourhood based biliteral filter
# Nadav Horesh
###########################################################
#
# Convert double buffer to short
# >> After a comparison between the paraller and the serial versions
#
cimport numpy as np
import numpy as np
import cython
from cython.parallel import prange
cdef extern from "math.h":
double exp(double x)
cdef extern from "math.h":
float expf(float x)
float fabsf(float x)
cdef extern:
int abs(int x)
############ T J UNMODIFIED CODE #############
cdef extern from "numpy/arrayobject.h":
ctypedef extern class numpy.flatiter [object PyArrayIterObject]:
cdef int nd_m1
cdef np.npy_intp index, size
cdef np.ndarray ao
cdef char *dataptr
# This isn't exposed to the Python API.
# So we can't use the same approach we used to define flatiter
ctypedef struct PyArrayNeighborhoodIterObject:
int nd_m1
np.npy_intp index, size
np.PyArrayObject *ao # note the change from np.ndarray
char *dataptr
object PyArray_NeighborhoodIterNew(flatiter it, np.npy_intp* bounds,
int mode, np.ndarray fill_value)
int PyArrayNeighborhoodIter_Next(PyArrayNeighborhoodIterObject *it)
int PyArrayNeighborhoodIter_Reset(PyArrayNeighborhoodIterObject *it)
object PyArray_IterNew(object arr)
void PyArray_ITER_NEXT(flatiter it)
np.npy_intp PyArray_SIZE(np.ndarray arr)
cdef enum:
NPY_NEIGHBORHOOD_ITER_ZERO_PADDING,
NPY_NEIGHBORHOOD_ITER_ONE_PADDING,
NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING,
NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING,
NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING
np.import_array()
############ END OF T J UNMODIFIED CODE #############
cdef int GAUSS_SAMP = 32
cdef int GAUSS_IDX_MAX = GAUSS_SAMP - 1
class Cbilateral_filter(object):
'''
A fully cythonic simple bilateral filtering.
The class provides the bilaterl filter function to be called by
generic_filter.
initialization parameters:
spat_sig: The sigma of the spatial Gaussian filter
inten_sig: The sigma of the gray-levels Gaussian filter
filter_size: (int) The size of the spatial convolution kernel. If
not set, it is set to ~ 4*spat_sig.
'''
def __init__(self, float spat_sig, float inten_sig, filter_size=None):
if filter_size is not None and filter_size >= 2:
self.xy_size = int(filter_size)
else:
self.xy_size = int(round(spat_sig*4))
# Make filter size odd
self.xy_size += 1-self.xy_size%2
x = np.arange(self.xy_size, dtype=np.float32)
x = (x-x.mean())**2
#xy_ker: Spatial convolution kernel
self.xy_ker = np.exp(-np.add.outer(x,x)/(2*spat_sig**2)).ravel()
self.xy_ker /= self.xy_ker.sum()
self.inten_sig = 2 * inten_sig**2
# self.index in the coordinate of the middle point
self.index = (self.xy_size+1) * (self.xy_size // 2)
## An initialization for LUT instead of a Gaussian function call
## (for the fc_filter method)
x = np.linspace(0,3.0, GAUSS_SAMP)
self.gauss_lut = np.exp(-x**2/2)
self.x_quant = 3*inten_sig / GAUSS_IDX_MAX
self.gauss_lut_float32 = np.exp(-x**2/2).astype(np.float32)
# @cython.boundscheck(False)
# @cython.wraparound(False)
@cython.cdivision(True)
def filterf(self, np.ndarray[dtype=np.float32_t, ndim=2] data not None):
# Define an iterator over the input array (image)
cdef arr_iter = PyArray_IterNew(<object>data)
# xsize, ysize: Input array dimensions
cdef unsigned int ysize = data.shape[0], xsize=data.shape[1]
# Define output array and a c pointer to iterate over it
cdef np.ndarray[np.float32_t, ndim=2] output = np.empty_like(data)
cdef float *out_ptr = <float *>output.data
# Get the alreasy initialized spatial and "z" Gaussian kernels
cdef np.ndarray[dtype=np.float32_t, ndim=1] kernel = self.xy_ker
cdef np.ndarray[dtype=np.float32_t, ndim=1] gauss_lut_arr =
self.gauss_lut_float32
# Iterators over the spatial kernel and input data
cdef float *pdata = <float *>data.data, *pker=<float *>kernel.data
cdef float *gauss_lut = <float *>gauss_lut_arr.data # C pointer to
iterate over "z"
# Misc temporary data for the convolution
cdef float sigma = self.inten_sig
cdef float weight_i, weight, result, centre, dat_i
cdef float x_quant = self.x_quant
# Misc untegers, loop indeces and spatial kernel size
cdef unsigned int i, j, k, centre_idx, exp_i
cdef int dim = self.xy_size
cdef unsigned int ker_lentgh = dim*dim
# Neighbourhood iteraor
cdef np.npy_intp* bounds_ptr = [-dim//2, dim//2, -dim//2, dim//2]
# Create the Python object and keep a reference to it
cdef object neigh_iter = PyArray_NeighborhoodIterNew(arr_iter,
bounds_ptr, NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING, None)
cdef PyArrayNeighborhoodIterObject *pneigh_iter = \
<PyArrayNeighborhoodIterObject *> neigh_iter
if not data.flags.c_contiguous:
raise ValueError('Input array must be C contiguous')
for i in range(ysize):
# for j in range(xsize):
for j in prange(xsize, nogil=True):
centre_idx = i*xsize+j
centre = pdata[centre_idx]
weight = 1.0E-30
result = 0.0
for k in range(ker_lentgh):
dat_i = (<float *>pneigh_iter.dataptr)[0]
exp_i = <int>(0.5 + fabsf((dat_i-centre) / x_quant))
weight_i = 0.0
if exp_i <= GAUSS_IDX_MAX:
weight_i = gauss_lut[exp_i] * pker[k]
weight += weight_i;
result += dat_i * weight_i
PyArrayNeighborhoodIter_Next(pneigh_iter)
out_ptr[centre_idx] = result/weight
PyArray_ITER_NEXT(arr_iter)
PyArrayNeighborhoodIter_Reset(pneigh_iter)
return output
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
