Introduction¶Image processing and analysis are generally seen as operations on two-dimensional arrays of values. There are however a number of fields where images of higher dimensionality must be analyzed. Good examples of these are medical imaging and biological imaging. Show
Table of Contents
Filter functions¶The functions described in this section all perform some type of spatial filtering of the input array: the elements in the output are some function of the values in the neighborhood of the corresponding input element. We refer to this neighborhood of elements as the filter kernel, which is often rectangular in shape but may also have an arbitrary footprint. Many of the functions described below allow you to define the footprint of the kernel, by passing a mask through the footprint parameter. For example a cross shaped kernel can be defined as follows: >>> footprint = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]]) >>> footprint array([[0, 1, 0], [1, 1, 1], [0, 1, 0]]) Usually the origin of the kernel is at the center calculated by dividing the dimensions of the kernel shape by two. For instance, the origin of a one-dimensional kernel of length three is at the second element. Take for example the correlation of a one-dimensional array with a filter of length 3 consisting of ones: >>> from scipy.ndimage import correlate1d >>> a = [0, 0, 0, 1, 0, 0, 0] >>> correlate1d(a, [1, 1, 1]) array([0, 0, 1, 1, 1, 0, 0]) Sometimes it is convenient to choose a different origin for the kernel. For this reason most functions support the origin parameter which gives the origin of the filter relative to its center. For example: >>> a = [0, 0, 0, 1, 0, 0, 0] >>> correlate1d(a, [1, 1, 1], origin = -1) array([0, 1, 1, 1, 0, 0, 0]) The effect is a shift of the result towards the left. This feature will not be needed very often, but it may be useful especially for filters that have an even size. A good example is the calculation of backward and forward differences: >>> a = [0, 0, 1, 1, 1, 0, 0] >>> correlate1d(a, [-1, 1]) # backward difference array([ 0, 0, 1, 0, 0, -1, 0]) >>> correlate1d(a, [-1, 1], origin = -1) # forward difference array([ 0, 1, 0, 0, -1, 0, 0]) We could also have calculated the forward difference as follows: >>> correlate1d(a, [0, -1, 1]) array([ 0, 1, 0, 0, -1, 0, 0]) However, using the origin parameter instead of a larger kernel is more efficient. For multidimensional kernels origin can be a number, in which case the origin is assumed to be equal along all axes, or a sequence giving the origin along each axis. Since the output elements are a function of elements in the neighborhood of the input elements, the borders of the array need to be dealt with appropriately by providing the values outside the borders. This is done by assuming that the arrays are extended beyond their boundaries according certain boundary conditions. In the functions described below, the boundary conditions can be selected using the mode parameter which must be a string with the name of the boundary condition. The following boundary conditions are currently supported:
The “constant” mode is special since it needs an additional parameter to specify the constant value that should be used. Note The easiest way to implement such boundary conditions would be to copy the data to a larger array and extend the data at the borders according to the boundary conditions. For large arrays and large filter kernels, this would be very memory consuming, and the functions described below therefore use a different approach that does not require allocating large temporary buffers. Correlation and convolution¶
Smoothing filters¶
Filters based on order statistics¶
Derivatives¶Derivative filters can be
constructed in several ways. The function
The Laplace filter is calculated by the sum of the second derivatives along all axes. Thus, different Laplace filters can be constructed using different second derivative functions. Therefore we provide a general function that takes a function argument to calculate the second derivative along a given direction.
The following two functions are
implemented using
The gradient magnitude is defined as the square root of the sum of the squares of the gradients in all directions. Similar to the generic Laplace function there is a
The
Generic filter functions¶To implement filter functions, generic functions can be used that accept a callable object that implements the filtering operation. The iteration over the input and output arrays is handled by these generic functions, along with such details as the implementation of the boundary
conditions. Only a callable object implementing a callback function that does the actual filtering work must be provided. The callback function can also be written in C and passed using a
These functions iterate over the lines or elements starting at the last axis, i.e. the last index changes the fastest. This order of iteration is guaranteed for the case that it is important to adapt the filter depending on spatial location. Here is an example of using a
class that implements the filter and keeps track of the current coordinates while iterating. It performs the same filter operation as described above for >>> a = np.arange(12).reshape(3,4) >>> >>> class fnc_class: ... def __init__(self, shape): ... # store the shape: ... self.shape = shape ... # initialize the coordinates: ... self.coordinates = [0] * len(shape) ... ... def filter(self, buffer): ... result = (buffer * np.array([1, 3])).sum() ... print self.coordinates ... # calculate the next coordinates: ... axes = range(len(self.shape)) ... axes.reverse() ... for jj in axes: ... if self.coordinates[jj] < self.shape[jj] - 1: ... self.coordinates[jj] += 1 ... break ... else: ... self.coordinates[jj] = 0 ... return result ... >>> fnc = fnc_class(shape = (3,4)) >>> generic_filter(a, fnc.filter, footprint = [[1, 0], [0, 1]]) [0, 0] [0, 1] [0, 2] [0, 3] [1, 0] [1, 1] [1, 2] [1, 3] [2, 0] [2, 1] [2, 2] [2, 3] array([[ 0, 3, 7, 11], [12, 15, 19, 23], [28, 31, 35, 39]]) For the >>> a = np.arange(12).reshape(3,4) >>> >>> class fnc1d_class: ... def __init__(self, shape, axis = -1): ... # store the filter axis: ... self.axis = axis ... # store the shape: ... self.shape = shape ... # initialize the coordinates: ... self.coordinates = [0] * len(shape) ... ... def filter(self, iline, oline): ... oline[...] = iline[:-2] + 2 * iline[1:-1] + 3 * iline[2:] ... print self.coordinates ... # calculate the next coordinates: ... axes = range(len(self.shape)) ... # skip the filter axis: ... del axes[self.axis] ... axes.reverse() ... for jj in axes: ... if self.coordinates[jj] < self.shape[jj] - 1: ... self.coordinates[jj] += 1 ... break ... else: ... self.coordinates[jj] = 0 ... >>> fnc = fnc1d_class(shape = (3,4)) >>> generic_filter1d(a, fnc.filter, 3) [0, 0] [1, 0] [2, 0] array([[ 3, 8, 14, 17], [27, 32, 38, 41], [51, 56, 62, 65]]) Fourier domain filters¶The functions described
in this section perform filtering operations in the Fourier domain. Thus, the input array of such a function should be compatible with an inverse Fourier transform function, such as the functions from the
Interpolation functions¶This section describes various interpolation functions that are based on B-spline theory. A good introduction to B-splines can be found in [1]. Spline pre-filters¶Interpolation using splines of an order larger than 1 requires a pre-filtering step. The interpolation functions described in section Interpolation functions apply pre-filtering by calling
Interpolation functions¶Following functions all employ spline interpolation to effect some type of geometric transformation of the input array. This requires a mapping of the output coordinates to the input coordinates, and therefore the possibility arises that input values outside the boundaries are needed. This problem is solved in the same way as described in Filter functions for the multidimensional filter functions. Therefore these functions all support a mode parameter that determines how the boundaries are handled, and a cval parameter that gives a constant value in case that the ‘constant’ mode is used.
Morphology¶Binary morphology¶
Most binary morphology functions can be expressed in terms of the basic operations erosion and dilation.
Here is an example of using >>> struct = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]]) >>> a = np.array([[1,0,0,0,0], [1,1,0,1,0], [0,0,1,1,0], [0,0,0,0,0]]) >>> a array([[1, 0, 0, 0, 0], [1, 1, 0, 1, 0], [0, 0, 1, 1, 0], [0, 0, 0, 0, 0]]) >>> from scipy.ndimage import binary_dilation >>> binary_dilation(np.zeros(a.shape), struct, -1, a, border_value=1) array([[ True, False, False, False, False], [ True, True, False, False, False], [False, False, False, False, False], [False, False, False, False, False]], dtype=bool) The
Other morphology operations can be defined in terms of erosion and d dilation. The following functions provide a few of these operations for convenience:
Grey-scale morphology¶Grey-scale morphology operations are the equivalents of binary morphology operations that operate on arrays with arbitrary values. Below we describe the grey-scale equivalents of erosion, dilation, opening and closing. These operations are implemented in a similar fashion as the filters described in Filter functions, and we refer to this section for the description of filter kernels and footprints, and the handling of array borders. The grey-scale morphology operations optionally take a structure parameter that gives the values of the structuring element. If this parameter is not given the structuring element is assumed to be flat with a value equal to zero. The shape of the structure can optionally be defined by the footprint parameter. If this parameter is not given, the structure is assumed to be rectangular, with sizes equal to the dimensions of the structure array, or by the size parameter if structure is not given. The size parameter is only used if both structure and footprint are not given, in which case the structuring element is assumed to be rectangular and flat with the dimensions given by size. The size parameter, if provided, must be a sequence of sizes or a single number in which case the size of the filter is assumed to be equal along each axis. The footprint parameter, if provided, must be an array that defines the shape of the kernel by its non-zero elements. Similar to binary erosion and dilation there are operations for grey-scale erosion and dilation:
Grey-scale opening and closing operations can be defined similar to their binary counterparts:
Distance transforms¶Distance transforms are used to calculate the minimum distance from each element of an object to the background. The following functions implement distance transforms for three different distance metrics: Euclidean, City Block, and Chessboard distances.
Segmentation and labeling¶Segmentation is the process of separating objects of interest from the background. The most simple approach is probably intensity thresholding, which is easily done with >>> a = np.array([[1,2,2,1,1,0], ... [0,2,3,1,2,0], ... [1,1,1,3,3,2], ... [1,1,1,1,2,1]]) >>> np.where(a > 1, 1, 0) array([[0, 1, 1, 0, 0, 0], [0, 1, 1, 0, 1, 0], [0, 0, 0, 1, 1, 1], [0, 0, 0, 0, 1, 0]]) The result is a binary image, in which the individual objects
still need to be identified and labeled. The function
There is a large number of other approaches for segmentation, for
instance from an estimation of the borders of the objects that can be obtained for instance by derivative filters. One such an approach is watershed segmentation. The function
Object measurements¶Given an array of labeled objects, the properties of the individual objects can be measured. The
The list of slices generated by >>> image = np.arange(4 * 6).reshape(4, 6) >>> mask = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]]) >>> labels = label(mask)[0] >>> slices = find_objects(labels) Then we can calculate the sum of the elements in the second object: >>> np.where(labels[slices[1]] == 2, image[slices[1]], 0).sum() 80 That is however not particularly efficient, and may also be more complicated for other types of measurements. Therefore a few measurements functions are defined that accept the array of object labels and the index of the object to be measured. For instance calculating the sum of the intensities can be done by: >>> from scipy.ndimage import sum as ndi_sum >>> ndi_sum(image, labels, 2) 80 For large arrays and small objects it is more efficient to call the measurement functions after slicing the array: >>> ndi_sum(image[slices[1]], labels[slices[1]], 2) 80 Alternatively, we can do the measurements for a number of labels with a single function call, returning a list of results. For instance, to measure the sum of the values of the background and the second object in our example we give a list of labels: >>> ndi_sum(image, labels, [0, 2]) array([178.0, 80.0]) The measurement functions described below all support the index parameter to indicate which object(s) should be measured. The default value of index is None. This indicates that all elements where the label is larger than zero should be treated as a single object and measured. Thus, in this case the labels array is treated as a mask defined by the elements that are larger than zero. If index is a number or a sequence of numbers it gives the labels of the objects that are measured. If index is a sequence, a list of the results is returned. Functions that return more than one result, return their result as a tuple if index is a single number, or as a tuple of lists, if index is a sequence.
Extending scipy.ndimage in C¶A few functions in An example of a function that supports callbacks is from scipy import ndimage def transform(output_coordinates, shift): input_coordinates = output_coordinates[0] - shift, output_coordinates[1] - shift return input_coordinates im = np.arange(12).reshape(4, 3).astype(np.float64) shift = 0.5 print(ndimage.geometric_transform(im, transform, extra_arguments=(shift,))) We can also implement the callback function with the following C code. /* example.c */ #include <Python.h> #include <numpy/npy_common.h> static int _transform(npy_intp *output_coordinates, double *input_coordinates, int output_rank, int input_rank, void *user_data) { npy_intp i; double shift = *(double *)user_data; for (i = 0; i < input_rank; i++) { input_coordinates[i] = output_coordinates[i] - shift; } return 1; } static char *transform_signature = "int (npy_intp *, double *, int, int, void *)"; static PyObject * py_get_transform(PyObject *obj, PyObject *args) { if (!PyArg_ParseTuple(args, "")) return NULL; return PyCapsule_New(_transform, transform_signature, NULL); } static PyMethodDef ExampleMethods[] = { {"get_transform", (PyCFunction)py_get_transform, METH_VARARGS, ""}, {NULL, NULL, 0, NULL} }; /* Initialize the module */ #if PY_VERSION_HEX >= 0x03000000 static struct PyModuleDef example = { PyModuleDef_HEAD_INIT, "example", NULL, -1, ExampleMethods, NULL, NULL, NULL, NULL }; PyMODINIT_FUNC PyInit_example(void) { return PyModule_Create(&example); } #else PyMODINIT_FUNC initexample(void) { Py_InitModule("example", ExampleMethods); } #endif More information on writing Python extension modules can be found here. If the C code is in the file from distutils.core import setup, Extension import numpy shift = Extension('example', ['example.c'], include_dirs=[numpy.get_include()] ) setup(name='example', ext_modules=[shift] ) and now running the script import ctypes import numpy as np from scipy import ndimage, LowLevelCallable from example import get_transform shift = 0.5 user_data = ctypes.c_double(shift) ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p) callback = LowLevelCallable(transform(), ptr) im = np.arange(12).reshape(4, 3).astype(np.float64) print(ndimage.geometric_transform(im, callback)) produces the same result as the original python script. In the C version The function
C callback functions for Below, we show alternative ways to write the code, using Cython, ctypes, or cffi instead of writing wrapper code in C. Numba Numba provides a way to write low-level functions easily in Python. We can write the above using Numba as: # example.py import numpy as np import ctypes from scipy import ndimage, LowLevelCallable from numba import cfunc, types, carray @cfunc(types.intc(types.CPointer(types.intp), types.CPointer(types.double), types.intc, types.intc, types.voidptr)) def transform(output_coordinates_ptr, input_coordinates_ptr, output_rank, input_rank, user_data): input_coordinates = carray(input_coordinates_ptr, (input_rank,)) output_coordinates = carray(output_coordinates_ptr, (output_rank,)) shift = carray(user_data, (1,), types.double)[0] for i in range(input_rank): input_coordinates[i] = output_coordinates[i] - shift return 1 shift = 0.5 # Then call the function user_data = ctypes.c_double(shift) ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p) callback = LowLevelCallable(transform.ctypes, ptr) im = np.arange(12).reshape(4, 3).astype(np.float64) print(ndimage.geometric_transform(im, callback)) Cython Functionally the same code as above can be written in Cython with somewhat less boilerplate as follows. # example.pyx from numpy cimport npy_intp as intp cdef api int transform(intp *output_coordinates, double *input_coordinates, int output_rank, int input_rank, void *user_data): cdef intp i cdef double shift = (<double *>user_data)[0] for i in range(input_rank): input_coordinates[i] = output_coordinates[i] - shift return 1 # script.py import ctypes import numpy as np from scipy import ndimage, LowLevelCallable import example shift = 0.5 user_data = ctypes.c_double(shift) ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p) callback = LowLevelCallable.from_cython(example, "transform", ptr) im = np.arange(12).reshape(4, 3).astype(np.float64) print(ndimage.geometric_transform(im, callback)) cffi With cffi, you can interface with a C function residing in a shared library (DLL). First, we need to write the shared library, which we do in C — this example is for Linux/OSX: /* example.c Needs to be compiled with "gcc -std=c99 -shared -fPIC -o example.so example.c" or similar */ #include <stdint.h> int _transform(intptr_t *output_coordinates, double *input_coordinates, int output_rank, int input_rank, void *user_data) { int i; double shift = *(double *)user_data; for (i = 0; i < input_rank; i++) { input_coordinates[i] = output_coordinates[i] - shift; } return 1; } The Python code calling the library is: import os import numpy as np from scipy import ndimage, LowLevelCallable import cffi # Construct the FFI object, and copypaste the function declaration ffi = cffi.FFI() ffi.cdef(""" int _transform(intptr_t *output_coordinates, double *input_coordinates, int output_rank, int input_rank, void *user_data); """) # Open library lib = ffi.dlopen(os.path.abspath("example.so")) # Do the function call user_data = ffi.new('double *', 0.5) callback = LowLevelCallable(lib._transform, user_data) im = np.arange(12).reshape(4, 3).astype(np.float64) print(ndimage.geometric_transform(im, callback)) You can find more information in the cffi documentation. ctypes With ctypes, the C code and the compilation of the so/DLL is as for cffi above. The Python code is different: # script.py import os import ctypes import numpy as np from scipy import ndimage, LowLevelCallable lib = ctypes.CDLL(os.path.abspath('example.so')) shift = 0.5 user_data = ctypes.c_double(shift) ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p) # Ctypes has no built-in intptr type, so override the signature # instead of trying to get it via ctypes callback = LowLevelCallable(lib._transform, ptr, "int _transform(intptr_t *, double *, int, int, void *)") # Perform the call im = np.arange(12).reshape(4, 3).astype(np.float64) print(ndimage.geometric_transform(im, callback)) You can find more information in the ctypes documentation. References¶
|