ENH: Use Cython to compute GLCM

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 + <int>(sin(angle) * distance + 0.5)
+                    col = c + <int>(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

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)):

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())
+          )