Merge pull request #2095 from ThomasWalter/cooc

Add uint16 images support for co-occurrence matrix

skimage/feature/_texture.pyx

@@ -7,13 +7,23 @@ cimport numpy as cnp
 from libc.math cimport sin, cos, abs
 from .._shared.interpolation cimport bilinear_interpolation, round
 from .._shared.transform cimport integrate
-
+import cython
 
 cdef extern from "numpy/npy_math.h":
     double NAN "NPY_NAN"
 
+ctypedef fused any_int:
+    cnp.uint8_t
+    cnp.uint16_t
+    cnp.uint32_t
+    cnp.uint64_t
+    cnp.int8_t
+    cnp.int16_t
+    cnp.int32_t
+    cnp.int64_t
 
-def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances,
+
+def _glcm_loop(any_int[:, ::1] image, double[:] distances,
                double[:] angles, Py_ssize_t levels,
                cnp.uint32_t[:, :, :, ::1] out):
     """Perform co-occurrence matrix accumulation.
@@ -21,15 +31,20 @@ def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances,
     Parameters
     ----------
     image : ndarray
-        Input image, which is converted to the uint8 data type.
+        Integer typed input image. Only positive valued images are supported. 
+        If type is other than uint8, the argument `levels` needs to be set.
     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],
+        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). This argument is required for
+        16-bit images or higher and is typically the maximum of the image. 
+        As the output matrix is at least `levels` x `levels`, it might
+        be preferable to use binning of the input image rather than 
+        large values for `levels`. 
     out : ndarray
         On input a 4D array of zeros, and on output it contains
         the results of the GLCM computation.
@@ -38,7 +53,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
+        any_int i, j
         cnp.float64_t angle, distance
 
     with nogil:

skimage/feature/tests/test_texture.py

@@ -15,6 +15,7 @@ class TestGLCM():
                                [0, 2, 2, 2],
                                [2, 2, 3, 3]], dtype=np.uint8)
 
+
     @test_parallel()
     def test_output_angles(self):
         result = greycomatrix(self.image, [1], [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4], 4)
@@ -50,6 +51,34 @@ class TestGLCM():
                              [0, 0, 2, 0]], dtype=np.uint32)
         np.testing.assert_array_equal(result[:, :, 0, 0], expected)
 
+    def test_error_raise_float(self):
+        for dtype in [np.float, np.double, np.float16, np.float32, np.float64]: 
+            np.testing.assert_raises(ValueError, greycomatrix, self.image.astype(dtype), [1], [np.pi], 4)
+
+    def test_error_raise_int_types(self):
+        for dtype in [np.int16, np.int32, np.int64, np.uint16, np.uint32, np.uint64]: 
+            np.testing.assert_raises(ValueError, greycomatrix, self.image.astype(dtype), [1], [np.pi])
+
+    def test_error_raise_negative(self):
+        np.testing.assert_raises(ValueError, greycomatrix, self.image.astype(np.int16) - 1, [1], [np.pi], 4)
+
+    def test_error_raise_levels_smaller_max(self):
+        np.testing.assert_raises(ValueError, greycomatrix, self.image - 1, [1], [np.pi], 3)
+        
+    def test_image_data_types(self):
+        for dtype in [np.uint16, np.uint32, np.uint64, np.int16, np.int32, np.int64]: 
+            img = self.image.astype(dtype)
+            result = greycomatrix(img, [1], [np.pi / 2], 4,
+                                  symmetric=True)
+            assert result.shape == (4, 4, 1, 1)
+            expected = np.array([[6, 0, 2, 0],
+                                 [0, 4, 2, 0],
+                                 [2, 2, 2, 2],
+                                 [0, 0, 2, 0]], dtype=np.uint32)
+            np.testing.assert_array_equal(result[:, :, 0, 0], expected)
+            
+        return
+    
     def test_output_distance(self):
         im = np.array([[0, 0, 0, 0],
                        [1, 0, 0, 1],

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.
 
@@ -20,18 +20,21 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False,
 
     Parameters
     ----------
-    image : array_like of uint8
-        Integer typed input image. The image will be cast to uint8, so
-        the maximum value must be less than 256.
+    image : array_like
+        Integer typed input image. Only positive valued images are supported. 
+        If type is other than uint8, the argument `levels` needs to be set.
     distances : array_like
         List of pixel pair distance offsets.
     angles : array_like
         List of pixel pair angles in radians.
     levels : int, optional
-        The input image should contain integers in [0, levels-1],
+        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). This argument is required for
+        16-bit images or higher and is typically the maximum of the image. 
+        As the output matrix is at least `levels` x `levels`, it might
+        be preferable to use binning of the input image rather than 
+        large values for `levels`. 
     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
@@ -50,7 +53,8 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False,
         `P[i,j,d,theta]` is the number of times that grey-level `j`
         occurs at a distance `d` and at an angle `theta` from
         grey-level `i`. If `normed` is `False`, the output is of
-        type uint32, otherwise it is float64.
+        type uint32, otherwise it is float64. The dimensions are:
+        levels x levels x number of distances x number of angles.
 
     References
     ----------
@@ -97,11 +101,28 @@ 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 np.issubdtype(image.dtype, np.float):
+        raise ValueError("Float images are not supported by greycomatrix. "
+                         "The image needs to be cast to an unsigned integer type.")
+    
+    # for image type > 8bit, levels must be set.
+    if image.dtype not in (np.uint8, np.int8) and levels is None:         
+        raise ValueError("The levels argument is required for data types other than uint8. "
+                         "The resulting matrix will be at least levels ** 2 in size.")
+                             
+    if image.dtype in (np.int8, np.int16, np.int32, np.int64) and np.any(image < 0):
+        raise ValueError("Negative valued images are not supported.")
+
+    if levels is None:    
+        levels = 256
+
+    if image_max >= levels:
+        raise ValueError("The image maximum needs to be smaller than `levels`.")
+
     distances = np.ascontiguousarray(distances, dtype=np.float64)
     angles = np.ascontiguousarray(angles, dtype=np.float64)