co-occurrence matrix can be calculated for 16 bit images

skimage/feature/_texture.pyx

@@ -13,7 +13,7 @@ cdef extern from "numpy/npy_math.h":
     double NAN "NPY_NAN"
 
 
-def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances,
+def _glcm_loop(cnp.uint16_t[:, ::1] image, double[:] distances,
                double[:] angles, Py_ssize_t levels,
                cnp.uint32_t[:, :, :, ::1] out):
     """Perform co-occurrence matrix accumulation.
@@ -29,7 +29,10 @@ def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances,
     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)
+        (typically 256 for an 8-bit image). Be aware that the co-occurrence
+        matrix for every angle and every distance is of size levels x levels. 
+        Choosing a too large level might result in exceedingly large matrix.
+        If you have 16 bit or 32 bit images, consider binning. 
     out : ndarray
         On input a 4D array of zeros, and on output it contains
         the results of the GLCM computation.
@@ -38,7 +41,7 @@ def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances,
 
     cdef:
         Py_ssize_t a_idx, d_idx, r, c, rows, cols, row, col
-        cnp.uint8_t i, j
+        cnp.uint16_t i, j
         cnp.float64_t angle, distance
 
     with nogil:

skimage/feature/texture.py

@@ -11,7 +11,7 @@ from ._texture import (_glcm_loop,
                        _multiblock_lbp)
 
 
-def greycomatrix(image, distances, angles, levels=256, symmetric=False,
+def greycomatrix(image, distances, angles, levels=None, symmetric=False,
                  normed=False):
     """Calculate the grey-level co-occurrence matrix.
 
@@ -30,8 +30,10 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False,
     levels : int, optional
         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). The maximum value is
-        256.
+        (typically 256 for an 8-bit image). Be aware that the co-occurrence
+        matrix for every angle and every distance is of size levels x levels. 
+        Choosing a too large level might result in exceedingly large matrix.
+        If you have 16 bit or 32 bit images, consider binning. 
     symmetric : bool, optional
         If True, the output matrix `P[:, :, d, theta]` is symmetric. This
         is accomplished by ignoring the order of value pairs, so both
@@ -97,11 +99,21 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False,
     assert_nD(distances, 1, 'distances')
     assert_nD(angles, 1, 'angles')
 
-    assert levels <= 256
     image = np.ascontiguousarray(image)
     assert image.min() >= 0
-    assert image.max() < levels
-    image = image.astype(np.uint8)
+
+    image_max = image.max()
+    if levels is None:
+        # if levels is not given, we assume that there are [0, image_max -1] levels to 
+        # be considered. There would be only zeros for all other rows and columns. 
+        levels = image_max + 1
+        
+    assert image_max < levels
+
+    # we cast to uint16 (because of fixed typing in cython) 
+    # this has no impact on the size of the co-occurrence matrix. 
+    image = image.astype(np.uint16)
+
     distances = np.ascontiguousarray(distances, dtype=np.float64)
     angles = np.ascontiguousarray(angles, dtype=np.float64)