diff --git a/skimage/feature/_greycomatrix.pyx b/skimage/feature/_greycomatrix.pyx new file mode 100644 index 00000000..b9ff7f23 --- /dev/null +++ b/skimage/feature/_greycomatrix.pyx @@ -0,0 +1,67 @@ +"""Cython implementation for computing a grey level co-occurance matrix +""" + +import numpy as np +cimport numpy as np +cimport cython + +cdef extern from "math.h": + double sin(double) + double cos(double) + +@cython.boundscheck(False) +def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2, + negative_indices=False, mode='c'] image, + np.ndarray[dtype=np.float64_t, ndim=1, + negative_indices=False, mode='c'] distances, + np.ndarray[dtype=np.float64_t, ndim=1, + negative_indices=False, mode='c'] angles, + int levels, + np.ndarray[dtype=np.uint32_t, ndim=4, + negative_indices=False, mode='c'] out + ): + """Perform co-occurnace matrix accumulation + + Parameters + ---------- + image : ndarray + Input image, which is converted to the uint8 data type. + distances : ndarray + List of pixel pair distance offsets. + angles : ndarray + List of pixel pair angles in radians. + levels : int + The input image should contain integers in [0, levels-1], + where levels indicate the number of grey-levels counted + (typically 256 for an 8-bit image) + out : ndarray + On input a 4D array of zeros, and on output it contains + the results of the GLCM computation. + + """ + cdef: + np.int32_t a_inx, d_idx + np.int32_t r, c, rows, cols, row, col + np.int32_t i, j + + rows = image.shape[0] + cols = image.shape[1] + + for a_idx, angle in enumerate(angles): + for d_idx, distance in enumerate(distances): + for r in range(rows): + for c in range(cols): + i = image[r, c] + + # compute the location of the offset pixel + row = r + (sin(angle) * distance + 0.5) + col = c + (cos(angle) * distance + 0.5); + + # make sure the offset is within bounds + if row >= 0 and row < rows and \ + col >= 0 and col < cols: + j = image[row, col] + + if i >= 0 and i < levels and \ + j >= 0 and j < levels: + out[i, j, d_idx, a_idx] += 1 diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 98631e28..b282dcff 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -6,6 +6,8 @@ properties to characterize image textures. import numpy as np import skimage.util +from _greycomatrix import _glcm_loop + def compute_glcm(image, distances, angles, levels=256, symmetric=False, normed=False): @@ -19,9 +21,9 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, image : ndarray Input image, which is converted to the uint8 data type. distances : array_like - List of histogram distance offsets. + List of pixel pair distance offsets. angles : array_like - List of histogram angles in radians. + List of pixel pair angles in radians. levels : int, optional The input image should contain integers in [0, levels-1], where levels indicate the number of grey-levels counted @@ -77,8 +79,8 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, assert image.ndim == 2 assert image.min() >= 0 assert image.max() < levels - distances = np.asarray(distances) - angles = np.asarray(angles) + distances = np.ascontiguousarray(distances, dtype=np.float64) + angles = np.ascontiguousarray(angles, dtype=np.float64) assert distances.ndim == 1 assert angles.ndim == 1 @@ -86,24 +88,8 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, out = np.zeros((levels, levels, len(distances), len(angles)), dtype=np.uint32) - for a_idx, angle in enumerate(angles): - for d_idx, distance in enumerate(distances): - for r in range(rows): - for c in range(cols): - i = image[r, c] - - # compute the location of the offset pixel - row = r + int(np.round(np.sin(angle) * distance)) - col = c + int(np.round(np.cos(angle) * distance)) - - # make sure the offset is within bounds - if row >= 0 and row < rows and \ - col >= 0 and col < cols: - j = image[row, col] - - if i >= 0 and i < levels and \ - j >= 0 and j < levels: - out[i, j, d_idx, a_idx] += 1 + # count co-occurances + _glcm_loop(image, distances, angles, levels, out) # make each GLMC symmetric if symmetric: @@ -111,7 +97,7 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, for a in range(len(angles)): out[:, :, d, a] += out[:, :, d, a].transpose() - # normalize each GLMC individually + # normalize each GLMC if normed: out = out.astype(np.float64) for d in range(len(distances)): diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py new file mode 100644 index 00000000..d2ab7d22 --- /dev/null +++ b/skimage/feature/setup.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python + +import os +from skimage._build import cython + +base_path = os.path.abspath(os.path.dirname(__file__)) + +def configuration(parent_package='', top_path=None): + from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs + + config = Configuration('feature', parent_package, top_path) + config.add_data_dir('tests') + + cython(['_greycomatrix.pyx'], working_path=base_path) + + config.add_extension('_greycomatrix', sources=['_greycomatrix.c'], + include_dirs=[get_numpy_include_dirs()]) + + return config + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(maintainer = 'Scikits-image Developers', + author = 'Scikits-image Developers', + maintainer_email = 'scikits-image@googlegroups.com', + description = 'Features', + url = 'https://github.com/scikits-image/scikits-image', + license = 'SciPy License (BSD Style)', + **(configuration(top_path='').todict()) + )