Add GLCM / Haralick attributes

bruges/attribute/__init__.py

@@ -21,3 +21,5 @@
 from .complex import quadrature
 
 from .horizon import *
+
+from .glcm import glcm2d

bruges/attribute/glcm.py

@@ -0,0 +1,158 @@
+"""
+GLCM (gray level co-occurrence matrices)
+:copyright: 2022 Software Underground 
+:license: Apache 2.0
+"""
+
+from bruges.util.util import patchify, unpatchify
+from skimage.feature.texture import greycomatrix, greycoprops
+from skimage.transform import resize, warp
+import skimage
+from scipy.ndimage import convolve
+import numpy as np
+from tqdm import tqdm
+
+def main():
+    pass
+
+def glcm2d(img, vmin=None, vmax=None, levels=8, kernel_size=5, distance=1.0, angle=0.0, feature='contrast', method='skimage'):
+    '''
+    Compute statistical properties from GLCM (gray level co-occurrence matrices) on seismic data.
+    https://subsurfwiki.org/wiki/Haralick_texture
+    Algorithm :
+        1. Read data from an image, seismic horizon, or grid
+        2. Compute GLCMs in a moving window
+        3. Compute statistics of GLCMs, and put in new grid
+        4. Result: various attributes, such as entropy, energy, contrast
+
+    Source:
+    GLCM
+        Haralick, R (1979). Statistical and structural approaches to texture. Proceedings of the IEEE 67, 786–804.
+    GLCM using skimage
+        Stéfan van der Walt, Johannes L. Schönberger, Juan Nunez-Iglesias, François Boulogne, Joshua D. Warner, Neil Yager, Emmanuelle Gouillart, Tony Yu and the scikit-image contributors. scikit-image: Image processing in Python. PeerJ 2:e453 (2014) https://doi.org/10.7717/peerj.453
+        https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_glcm.html?highlight=glcm
+    GLCM on seismic data
+        Marcilio Castro de Matos, Malleswar (Moe) Yenugu, Sipuikinene Miguel Angelo, and Kurt J. Marfurt, (2011), "Integrated seismic texture segmentation and cluster analysis applied to channel delineation and chert reservoir characterization," GEOPHYSICS 76: P11-P21.https://doi.org/10.1190/geo2010-0150.1
+        https://library.seg.org/doi/abs/10.1190/geo2010-0150.1?journalCode=gpysa7#
+    Fast GLCM 
+        tzm030329, Fast GLCM, (2018), GitHub repository, https://github.com/tzm030329/GLCM
+
+    Args:
+    ----------
+    img: 2darray, shape=(h,w), dtype=np.uint8 / float
+        input image
+    vmin: float
+        minimum value of input image
+    vmax: float
+        maximum value of input image
+    levels: int
+        number of grey-levels of GLCM
+    kernel_size: int
+        Patch size to calculate GLCM around the target pixel
+    distance: float
+        pixel distance offsets (1.0, 2.0, and etc.)
+    angle: float
+        pixel angles (0.0, 30.0, 45.0, 90.0, and etc.)
+    feature: string
+        statistical feature from GLCM ('contrast','dissimilarity','homogeneity','ASM','energy')
+    method: string
+        method to use : ('skimage', 'fast')
+            'skimage': default,
+            'fast': Fastest GLCM calculation using opencv warp affine and 2D Convolution
+
+    Returns:
+    -------
+    result: 2darray
+        2D Grey-level co-occurrence matrix statistical feature array the same shape as the input array.
+
+    '''
+    if (vmin == None) or (vmax == None):
+        mi, ma = img.min(), img.max()
+    else:
+        mi, ma = vmin, vmax
+    ks = kernel_size
+    h,w = img.shape
+
+    # digitize
+    bins = np.linspace(mi, ma+1, levels+1)
+    gl1 = np.digitize(img, bins) - 1
+
+
+    if method == 'skimage':        
+        patches = patchify(gl1, (kernel_size, kernel_size), step=1)
+        glcmstat = np.zeros_like(patches[:,:,0,0])
+        with tqdm(total=patches.shape[0]*patches.shape[1], desc="GLCM statistic calculation") as pbar:
+            for i in range(patches.shape[0]):
+                for j in range(patches.shape[1]):
+                    glcm2 = greycomatrix(patches[i,j,:,:], distances=[int(distance)], angles=[int(angle)], levels=levels,
+                                            #symmetric=True, 
+                                            normed=True
+                                            )
+                    result2 = greycoprops(glcm2, feature)
+                    glcmstat[i,j] = result2[0][0]
+                    pbar.update(1)
+        pbar.close()        
+        result = resize(glcmstat,(gl1.shape[0], gl1.shape[1]))
+
+    if method == 'fast':
+        # make shifted image
+        dx = distance*np.cos(np.deg2rad(angle))
+        dy = distance*np.sin(np.deg2rad(-angle))
+        #mat = np.array([[1.0,0.0,-dx], [0.0,1.0,-dy]], dtype=np.float32)
+        mat2 = np.array([[1.0,0.0,-dx], [0.0,1.0,-dy], [0.0,0.0,1.0]], dtype=np.float32)
+        gl2 = warp(gl1.astype(float), mat2, mode='edge')
+        # make glcm
+        glcm = np.zeros((levels, levels, h, w), dtype=np.uint8)
+        for i in range(levels):
+            for j in range(levels):
+                mask = ((gl1==i) & (gl2==j))
+                glcm[i,j, mask] = 1
+
+        kernel = np.ones((ks, ks), dtype=np.uint8)
+        for i in range(levels):
+            for j in range(levels):
+                glcm[i,j] = convolve(glcm[i,j], kernel, mode='nearest')
+
+        glcm = glcm.astype(np.float32)
+
+        #select feature
+        if feature == 'contrast':
+            result = np.zeros((h,w), dtype=np.float32)
+            for i in range(levels):
+                for j in range(levels):
+                    result += glcm[i,j] * (i-j)**2 
+
+        if feature == 'dissimilarity':
+            result = np.zeros((h,w), dtype=np.float32)
+            for i in range(levels):
+                for j in range(levels):
+                    result += glcm[i,j] * np.abs(i-j)
+
+        if feature == 'homogeneity':    
+            result = np.zeros((h,w), dtype=np.float32)
+            for i in range(levels):
+                for j in range(levels):
+                    result += glcm[i,j] / (1.+(i-j)**2)
+
+        if feature == 'ASM':
+            result = np.zeros((h,w), dtype=np.float32)
+            for i in range(levels):
+                for j in range(levels):
+                    result  += glcm[i,j]**2
+
+        if feature == 'energy':
+            result = np.zeros((h,w), dtype=np.float32)
+            for i in range(levels):
+                for j in range(levels):
+                    result  += glcm[i,j]**2
+            result = np.sqrt(result)
+    #normalizing        
+    result =  result/(result.max()/1.0)        
+    return result
+
+
+if __name__ == '__main__':
+    main()
+
+    img = skimage.data.camera()
+    result = glcm2d(data, levels=8, kernel_size=3, distance=1.0, angle=0.0, feature='dissimilarity', method='fast')

bruges/util/util.py

@@ -11,6 +11,8 @@
 import scipy.signal
 import numpy as np
 
+from skimage.util.shape import view_as_windows
+from typing import Tuple, Union, cast
 
 def deprecated(instructions):
     """
@@ -309,3 +311,173 @@ def power(start, stop, num):
     x = np.linspace(0, 8, num)
     y = 1 - 2**-x
     return min(start, stop) + abs(stop-start) * y
+
+
+Imsize = Union[Tuple[int, int], Tuple[int, int, int]]
+
+def patchify(image: np.ndarray, patch_size: Imsize, step: int = 1) -> np.ndarray:
+    """
+    Split a 2D or 3D image into small patches given the patch size.
+    Parameters
+    ----------
+    image: the image to be split. It can be 2d (m, n) or 3d (k, m, n)
+    patch_size: the size of a single patch
+    step: the step size between patches
+    Examples
+    --------
+    >>> image = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
+    >>> patches = patchify(image, (2, 2), step=1)  # split image into 2*3 small 2*2 patches.
+    >>> assert patches.shape == (2, 3, 2, 2)
+    >>> reconstructed_image = unpatchify(patches, image.shape)
+    >>> assert (reconstructed_image == image).all()
+    """
+    return view_as_windows(image, patch_size, step)
+
+
+def unpatchify(patches: np.ndarray, imsize: Imsize) -> np.ndarray:
+    """
+    Merge patches into the orignal image
+    Parameters
+    ----------
+    patches: the patches to merge. It can be patches for a 2d image (k, l, m, n)
+             or 3d volume (i, j, k, l, m, n)
+    imsize: the size of the original image or volume
+    Examples
+    --------
+    >>> image = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
+    >>> patches = patchify(image, (2, 2), step=1)  # split image into 2*3 small 2*2 patches.
+    >>> assert patches.shape == (2, 3, 2, 2)
+    >>> reconstructed_image = unpatchify(patches, image.shape)
+    >>> assert (reconstructed_image == image).all()
+    """
+
+    assert len(patches.shape) / 2 == len(
+        imsize
+    ), "The patches dimension is not equal to the original image size"
+
+    if len(patches.shape) == 4:
+        return _unpatchify2d(patches, cast(Tuple[int, int], imsize))
+    elif len(patches.shape) == 6:
+        return _unpatchify3d(patches, cast(Tuple[int, int, int], imsize))
+    else:
+        raise NotImplementedError(
+            "Unpatchify only supports a matrix of 2D patches (k, l, m, n)"
+            f"or 3D volumes (i, j, k, l, m, n), but got: {patches.shape}"
+        )
+
+
+def _unpatchify2d(  # pylint: disable=too-many-locals
+    patches: np.ndarray, imsize: Tuple[int, int]
+) -> np.ndarray:
+
+    assert len(patches.shape) == 4
+
+    i_h, i_w = imsize
+    image = np.zeros(imsize, dtype=patches.dtype)
+
+    n_h, n_w, p_h, p_w = patches.shape
+
+    s_w = 0 if n_w <= 1 else (i_w - p_w) / (n_w - 1)
+    s_h = 0 if n_h <= 1 else (i_h - p_h) / (n_h - 1)
+
+    # The step size should be same for all patches, otherwise the patches are unable
+    # to reconstruct into a image
+    if int(s_w) != s_w:
+        raise NonUniformStepSizeError(i_w, n_w, p_w, s_w)
+    if int(s_h) != s_h:
+        raise NonUniformStepSizeError(i_h, n_h, p_h, s_h)
+    s_w = int(s_w)
+    s_h = int(s_h)
+
+    i, j = 0, 0
+
+    while True:
+        i_o, j_o = i * s_h, j * s_w
+
+        image[i_o : i_o + p_h, j_o : j_o + p_w] = patches[i, j]
+
+        if j < n_w - 1:
+            j = min((j_o + p_w) // s_w, n_w - 1)
+        elif i < n_h - 1 and j >= n_w - 1:
+            # Go to next row
+            i = min((i_o + p_h) // s_h, n_h - 1)
+            j = 0
+        elif i >= n_h - 1 and j >= n_w - 1:
+            # Finished
+            break
+        else:
+            raise RuntimeError("Unreachable")
+
+    return image
+
+
+def _unpatchify3d(  # pylint: disable=too-many-locals
+    patches: np.ndarray, imsize: Tuple[int, int, int]
+) -> np.ndarray:
+
+    assert len(patches.shape) == 6
+
+    i_h, i_w, i_c = imsize
+    image = np.zeros(imsize, dtype=patches.dtype)
+
+    n_h, n_w, n_c, p_h, p_w, p_c = patches.shape
+
+    s_w = 0 if n_w <= 1 else (i_w - p_w) / (n_w - 1)
+    s_h = 0 if n_h <= 1 else (i_h - p_h) / (n_h - 1)
+    s_c = 0 if n_c <= 1 else (i_c - p_c) / (n_c - 1)
+
+    # The step size should be same for all patches, otherwise the patches are unable
+    # to reconstruct into a image
+    if int(s_w) != s_w:
+        raise NonUniformStepSizeError(i_w, n_w, p_w, s_w)
+    if int(s_h) != s_h:
+        raise NonUniformStepSizeError(i_h, n_h, p_h, s_h)
+    if int(s_c) != s_c:
+        raise NonUniformStepSizeError(i_c, n_c, p_c, s_c)
+
+    s_w = int(s_w)
+    s_h = int(s_h)
+    s_c = int(s_c)
+
+    i, j, k = 0, 0, 0
+
+    while True:
+
+        i_o, j_o, k_o = i * s_h, j * s_w, k * s_c
+
+        image[i_o : i_o + p_h, j_o : j_o + p_w, k_o : k_o + p_c] = patches[i, j, k]
+
+        if k < n_c - 1:
+            k = min((k_o + p_c) // s_c, n_c - 1)
+        elif j < n_w - 1 and k >= n_c - 1:
+            j = min((j_o + p_w) // s_w, n_w - 1)
+            k = 0
+        elif i < n_h - 1 and j >= n_w - 1 and k >= n_c - 1:
+            i = min((i_o + p_h) // s_h, n_h - 1)
+            j = 0
+            k = 0
+        elif i >= n_h - 1 and j >= n_w - 1 and k >= n_c - 1:
+            # Finished
+            break
+        else:
+            raise RuntimeError("Unreachable")
+
+    return image
+
+
+class NonUniformStepSizeError(RuntimeError):
+    def __init__(
+        self, imsize: int, n_patches: int, patch_size: int, step_size: float
+    ) -> None:
+        super().__init__(imsize, n_patches, patch_size, step_size)
+        self.n_patches = n_patches
+        self.patch_size = patch_size
+        self.imsize = imsize
+        self.step_size = step_size
+
+    def __repr__(self) -> str:
+        return f"Unpatchify only supports reconstructing image with a uniform step size for all patches. \
+However, reconstructing {self.n_patches} x {self.patch_size}px patches to an {self.imsize} image requires {self.step_size} as step size, which is not an integer."
+
+    def __str__(self) -> str:
+        return self.__repr__()

requirements.txt

@@ -1 +1,4 @@
-.
+scikit-image>=0.16.2
+tqdm>=4.47.0
+numpy>=1.19.5
+scipy>=1.5.0