diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index 57d8ec81..fb4189c8 100644 --- a/skimage/feature/_texture.pyx +++ b/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: diff --git a/skimage/feature/tests/test_texture.py b/skimage/feature/tests/test_texture.py index ae54a6bb..bfa377e6 100644 --- a/skimage/feature/tests/test_texture.py +++ b/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], diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py index 4be6c749..4de65c76 100644 --- a/skimage/feature/texture.py +++ b/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)