Scaling laws for Haralick texture features of linear gradients

The Gray Level Co-occurrence Matrix

A grayscale image is a two-dimensional matrix $I(x,y)$ that stores gray-level intensities (see Fig. 1A). The bit depth of an image determines the number $N_g$ of gray levels. For instance, an 8-bit image has $N_g=2^8=256$ gray levels. By convention, a gray level of zero, $I(x,y)=0$, represents black, while $I(x,y)=N_g-1$ corresponds to white. Intermediate intensities represent various shades of gray. Figure 1A shows the upper left corner reference frame attached to an image with x-direction pointing horizontally to the right and the y-direction vertically downward. Each square in Fig. 1A represents an image pixel. Arrows from the central highlighted pixel indicate the offset vectors $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$ to its neighbors. The increment $\mathrm{\Delta }x$ represents the image row offset and $\mathrm{\Delta }y$ represents the image column offset.

Figure 1B illustrates a rectangular $N_x(=5)\times N_y(=4)$ image with a 2-bit depth ( $N_g\in \{0,1,2,3\}$). For the same image, as shown in Fig. 1B, each displacement vector $d$ defines a corresponding GLCM. For instance, a unit displacement along horizontal direction $d=(\mathrm{\Delta }x=1,\mathrm{\Delta }y=0)$ produces Fig. 1C. Indeed, there are two pairs of pixels with the starting point gray level $i=0$ and endpoint intensity level $j=1$ separated by one pixel displacement along the horizontal direction. The array coordinates (1,4)-(1,5) and (2,2)-(2,3) are marked with elliptical shaded area and connected by the two horizontal lines extending from panel B image to the corresponding GLCM entry $P(0,1)=2$ in panel C. Similarly, there is only one pair of pixels in the Fig. 1B image with the starting point gray level $i=1$ and endpoint intensity level $j=2$ separated by one pixel displacement along the vertical direction. The array coordinates (3,4)-(4,4) are marked with rectangular shaded area and connected horizontally by a line extending from panel B image to the corresponding GLCM entry $P(1,2)=1$ in panel D. The unnormalized GLCM counts the number of occurrences of the (reference) gray level $i$ at a distance specified by the displacement vector $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$ from the (target) gray level $j$ (Haralick, Shanmugam & Dinstein, 1973):

(1) $$P_d(i,j)=\mathrm{\#}\{((x_i,y_i),(x_j,y_j)):I(x_i,y_i)=i\mathrm{\&}I(x_j,y_j)=j\},$$where # denotes the number of elements in the set, the coordinates of the reference gray level $i$ are $(x_i,y_i)$, and the coordinates of the neighbor (target) pixel with gray level $j$ are $(x_j=x_i+\mathrm{\Delta }x,y_j=y_i+\mathrm{\Delta }y)$.

In the GLCM Eq. (1), the first index $i$ represents the intensity of the reference point, or the starting point of the displacement vector $d$, while the second index $j$ corresponds to the intensity of the endpoint of the displacement vector. For instance, an offset $d=(\mathrm{\Delta }x=1,\mathrm{\Delta }y=0)$ indicates that the row index (in the $y$-direction) remains unchanged since $\mathrm{\Delta }y=0$, and the column index (in the $x$- or horizontal direction across the image) increases by one unit $(\mathrm{\Delta }x=1)$.

For simplicity, Fig. 1 only counts the pairs for gray levels one pixel apart along the horizontal (Fig. 1C) and vertical (Fig. 1D) directions, respectively. For example, only one pair of gray level intensities 1-0 is counted between the spatial coordinates (2,3) and (2,4) in Fig. 1B, which is shown as $P(1,0)=1$ in Fig. 1C GLCM. As a result, the Fig. 1B GLCM is not symmetric. In the original definition of the GLCM provided by Haralick (Haralick, Shanmugam & Dinstein, 1973), symmetry allows both $P(1,2)$ and $P(2,1)$ pairings to be counted as instances where the pixel value 1 is separated by the distance vector $d$ from the pixel value 2. Mathematically this is achieved by adding to the GLCMs in Figs. 1C and 1D their corresponding transposed arrays. In line with Haralick’s definition, our implementation and all the results presented in this study used a symmetric GLCM matrix definition.

The number of possible pairs in the image typically normalizes the GLCM. For instance, in an $N_x\times N_y$ image, there are $R_x=(N_x-1)N_y$ horizontal pairs and $R_y=(N_y-1)N_x$ vertical pairs. In the example depicted in Fig. 1, since the image has $4\times 5$ pixels, the GLCM normalization factors are $R_x=16$ and $R_y=15$. The corresponding normalized GLCM values in Fig. 1C are, for example, $p_d(0,2)=P(0,2)/R_x=2/16$, and for Fig. 1D, they are $p_d(1,2)=P_d(1,2)/R_y=1/15$. The unnormalized GLCM is indicated with capital letters such as $P_d(i,j)$, while its normalized version is denoted as $p_d(i,j)$:

(2) $$p_d(i,j)=\frac{P_d(i,j)}{{\sum }_{i=0}^{N_g-1}{\sum }_{j=0}^{N_g-1}{P}_d(i,j)}.$$

The normalized GLCM indicates the likelihood of finding gray level $j$ at a displacement $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$ from the current location of the reference pixel with gray level $i$ in an image. It adheres to the normalization condition ${\sum }_{i=0}^{N_g-1}{\sum }_{j=0}^{N_g-1}{p}_d(i,j)=1$. More than half of the original 14 Haralick features rely on an additional step that involves computing marginal probability distributions from $p_d(i,j)$.

The GLCM is a natural measure of image gradients, quantifying the change in light intensity from the reference intensity $i$ to the target intensity $j$ along the displacement vector $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$. Since Haralick features are scalar measures defined by the two-point histogram represented by the GLCM, they also inherently measure light intensity gradients present in images.

Marginal distributions associated with the GLCM

Only three of the original Haralick features (Haralick, Shanmugam & Dinstein, 1973; Haralick, 1979) use the normalized GLCM $p_d(i,j)$ as defined in Eq. (2). All the other use one of the four marginal probability distributions derived from $p_d(i,j)$. To simplify the notation, one dropped the subscript $d$ from the normalized GLCM $p_d(i,j)$. The $x$-direction marginal probability distribution can be obtained by summing along the rows of the GLCM $p(i,j)$:

$$p_x(i)=\sum _{j=0}^{N_g-1}p(i,j),$$as shown in Fig. 2A. For example, $p_x(0)$ is the sum of all row elements with an intensity $i=0$ at the reference point (see Fig. 1A), regardless of the intensity of its endpoint determined by the displacement vector. Therefore, $p_x(i)$ gives the probability of finding gray level $i$ in the image. The mean and variance of the GLCM along the marginal distribution $p_x(i)$ are ${\mu }_x={\sum }_{i=0}^{N_g-1}i{p}_x(i),$ and ${\sigma }_x^2={\sum }_{i=0}^{N_g-1}{(i-{\mu }_x)}^2{p}_x(i).$

The $y$-direction marginal probability distribution $p_y(i)$ can be obtained by summing the columns of the GLCM $p(i,j)$:

$$p_y(j)=\sum _{i=0}^{N_g-1}p(i,j).$$

For example, $p_y(0)$ is the sum of all column elements with an endpoint intensity $j=0$, regardless of the intensity of the reference (starting) point. These marginal probabilities are illustrated in Fig. 2, along the horizontal dashed lines representing the GLCM line summation for $p_x$ and along the vertical dashed lines representing the GLCM column summation of $p(i,j)$ to obtain $p_y$, respectively.

The marginal distribution of gray level differences $k=i-j$ between the reference pixel intensity $i$ and the endpoint intensity $j$ determined by the displacement vector $d$ is:

(3) $$p_{x-y}(k)=\sum _{i=0}^{N_g-1}\sum _{j=0}^{N_g-1}{\delta }_{|i-j|,k}p(i,j),$$where ${\delta }_{m,n}$ is Kronecker’s symbol. For example, $p_{x-y}(0)$ represents the sum of all primary diagonal elements of the GLCM, as these elements exhibit no gray level differences between the reference point and the endpoint of the vector $d$, as illustrated in Fig. 2B. Similarly, the sum of the elements along the first line parallel to and above the primary diagonal reflects a gray level difference of $k=+1$ units between the reference gray level $i$ and the endpoint gray level $j$ of the GLCM, which defines $p_{x-y}(1)$. The sum $p(0,1)+p(1,2)+p(2,3)$ of GLCM entries along the first line parallel and above the primary diagonal in Fig. 2B correspond to the fraction of $p_{x-y}(1)$ with $j-i=+1$. The sum $p(1,0)+p(2,1)+p(3,2)$ of GLCM entries along the first line parallel and below the primary diagonal in Fig. 2B correspond to the fraction of $p_{x-y}(1)$ with $j-i=-1$. Since the definition of gray level differences marginal distribution $p_{x-y}(1)$ in Eq. (3) counts absolute differences $k=|i-j|$, the two partial sums must also be added (see the $\oplus$ symbol) to produce $p_{x-y}(1)$.

The marginal distribution of gray level sums $k=i+j$ between the reference pixel intensity $i$ and the endpoint neighbor intensity $j$ is:

(4) $$p_{x+y}(k)=\sum _{i=0}^{N_g-1}\sum _{j=0}^{N_g-1}{\delta }_{i+j,k}p(i,j).$$

To prevent overcrowding in Fig. 2B, we only showed the $p_{x+y}(3)$, which signifies the sum of the GLCM elements along its secondary diagonal with $i+j=3$, i.e., $p(3,0)+(2,1)+p(1,2)+p(0,3)$. Other values for $p_{x+y}(s)$ correspond to summation along lines parallel to the secondary diagonal in Fig. 2B.

Synthetic gradient images

While the GLCM method described in “Methods” applies to any image, this study specifically focuses on computer-generated (synthetic) images with one-dimensional vertical gradients. This focus is motivated by the fact that image gradients are highly invariant across images (Long & Purves, 2003; Tward, 2021; Dresp-Langley & Reeves, 2024).

Image gradients have long been used as statistical (or spectral) priors for estimating image features (Gong & Sbalzarini, 2014, 2016). A gradient image $G(x,y)$, is derived from the first-order spatial differences of the original image, $I(x,y)$, such that $G(x,y)=(I(x-1,y)-I(x,y),I(x,y-1)-I(x,y))$ (McCann & Pollard, 2008; Sevcenco & Agathoklis, 2021). The gradient image retains the same dimensions as the original but stores the x- and y-direction gradient values at each pixel.

Gradient spectral priors have been extensively applied in various image processing tasks, including denoising and deblurring (Chen, Yang & Wu, 2010), image restoration (Cho et al., 2012), range compression (Fattal, Lischinski & Werman, 2002), shadow removal (Finlayson, Hordley & Drew, 2002), and image compositing (Levin et al., 2004; Perez, Gangnet & Blake, 2003). Notably, deblurring in the gradient domain is often more computationally efficient than operating on raw pixel values (Cho & Lee, 2009; Shan, Jia & Agarwala, 2008; Wang & Cheng, 2016).

Traditionally, images are decomposed into 2D orthogonal gradient maps assuming that x- and y-direction gradients are statistically independent. One of the first studies to explore potential correlations between these gradient distributions in natural scene images found “weakly negatively correlated in the training dataset (from edges in the images)” (Gong & Sbalzarini, 2016). Consistent with these findings, recent algorithms for image denoising and deblurring (Zheng et al., 2022; Zhangying et al., 2024), range compression (Yan, Sun & Davis, 2024), or pattern classification (Wang et al., 2025) continue to treat orthogonal gradients as independent and their spectral priors as uncorrelated. Based on this well-supported assumption, our study focuses exclusively on a vertical gradient for calculating Haralick texture features.

Figure 3A shows a $b=3$-bit depth grayscale image with dimensions $N_x\times N_y$, featuring a vertical, linearly increasing, periodic intensity gradient of $\mathrm{\nabla }=1$ gray level per pixel. The array $I(x,y)$ that represents the image is given by $I(x_i,y_i)=y_i\mathrm{\nabla }$ where $y_i=\{0,2,...,N_y-1\}$. Since image intensities do not depend on the $x_i=\{0,2,...,N_x-1\}$ matrix index, the image appears as horizontal stripes with linearly increasing intensity (Fig. 3A). Furthermore, the vertical gradients are periodic, i.e., the intensity pattern repeats after reaching the maximum number of gray levels $N_g=2^b$. In other words, the vertical coordinate $y_i$ and pixel intensity are connected through $I(x_i,y_i)=mod(y_i,N_g)\mathrm{\nabla }$. The modulo (“mod”) operation along the vertical spatial indices $y_i$ ensures the gradient repeats periodically after $N_g$ pixels. In Fig. 3A example, the gray levels increase linearly from zero to $N_g-1$ with a step of $\mathrm{\nabla }=1$ gray level per pixel. The arrow next to the gradient in Fig. 3A indicates the gradient’s direction. Similarly, Fig. 3E shows a synthetically generated image with a vertical, linearly increasing, and periodic gradient $\mathrm{\nabla }=2$ gray levels per pixel. The grayscale images from Fig. 3A and Fig. 3E are numerically represented in Fig. 3B and Fig. 3F, respectively. The horizontal arrows between panels A and B indicate that the constant intensity line of pixels is represented numerically by the corresponding integer values with black mapped to 0. Following the procedure described above, we generated square synthetic images of $1024\times 1024$ pixels containing periodic linear gradient patterns, as illustrated in Fig. 3. Our analysis focuses on three key variables:

• (1)

  The number of gray levels in the image ( $N_g$),

• (2)

  The intensity of the image gradients ( $\mathrm{\nabla }$) in gray levels per pixel, and

• (3)

  The displacement vector ( $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$) in pixels, which determines the GLCM matrix used to compute the Haralick features.

We created images with a bit depth (b) ranging from 4 to 8, corresponding to $N_g=2^b\in \{16,32,64,128,256\}$. These values represent a broad and realistic range for evaluating how Haralick features depend on $N_g$ (see Figs. 4 and 5). For each bit depth we generated synthetic images with gradient intensities ( $\mathrm{\nabla }$) ranging from 1 to 8. However, to reduce visual clutter, only odd $\mathrm{\nabla }$ values are displayed in Figs. 4 and 5. Finally, for each combination of bit depth ( $b$) and gradient intensity $\mathrm{\nabla }$, we computed GLCMs for vertical displacement vectors $|d|=1,\dots ,8$.

Two-dimensional gradient maps

To count the pairs of pixels with a starting gray level $i$ and an endpoint gray level $j$ separated by a distance $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$ pixels, one can create a two-dimensional (2D) $N_y\times N_y$ gradient map such that its $(y_i,y_j)=(y_i,y_{j=i+|d|})$ entry is 1 if $I(x_i,y_{j=i+|d|})-I(x_i,y_i)=d\cdot \mathrm{\nabla }$ and zero otherwise as shown in Fig. 3C. Here, $\cdot$ is the dot product and ensures that one considers the relative orientation of the gradient $\mathrm{\nabla }$ to the displacement vector $d$. From Figs. 3B and 3F, one notices that the gray level intensity at spatial coordinate $y_i$ is always $y_i=i\mathrm{\nabla }$ with $i=0,\dots ,N_y-1$. The pixel intensity at a vertical coordinate $y_j$, which is a distance $|d|$ from $y_i$, is $y_j=y_i+d\cdot \mathrm{\nabla }=(i+|d|)\mathrm{\nabla }$. As a result, the 2D maps in Figs. 3C and 3G are one-to-one correspondences between pixel location $y_i$ and its corresponding gray level intensity $i\mathrm{\nabla }$. One notices in Fig. 3C with $\mathrm{\nabla }=1$ gray level per pixel and Fig. 3G with $\mathrm{\nabla }=2$ gray level per pixel that the vertical and horizontal distance between any non-zero entries of the 2D gradient map is $\mathrm{\nabla }$. These $\mathrm{\nabla }$ displacements are marked in Fig. 3C and Fig. 3G, respectively. Additionally, one can observe from the 2D gradient maps in Figs. 3C and 3G that all nonzero entries are aligned with the primary diagonal of the 2D gradient map at a distance of $d\cdot \mathrm{\nabla }$ from it. The distance of the gradient pattern from the primary diagonal of the 2D gradient maps is determined by the first gray level intensity, i.e., $i=0$, which is always paired with the gray label $j=d\cdot \mathrm{\nabla }$ for any displacement vector $d$ and gradient intensity $\mathrm{\nabla }$. Finally, all nonzero entries $(y_i,y_j)$ in the 2D gradient maps shown in Figs. 3C and 3G obey the condition $k=|i-j|=d\cdot \mathrm{\nabla }$ shown with dashed line parallel to the principal axis diagonal. The principal diagonal elements are always zero because they correspond to a uniform image with no intensity changes from pixel to pixel.

Wrap around the 2D gradient map to get the GLCM

While illuminating, representing a periodic linear gradient of length $N_y$ using a sparse $N_y\times N_y$ 2D gradient maps, as shown in Figs. 3C and 3G, is not efficient. As a result, the GLCM removes the extra spatial information about pixel coordinates $(y_i,y_j)$ retained by the 2D gradient map and only counts the co-occurrence of gray level intensities $i$ and $j$ at a relative distance $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$, as illustrated in Figs. 3D and 3H. Consequently, for a specific displacement vector $d$, the GLCM is an $N_g\times N_g$ matrix that solely counts the co-occurrence of gray levels $i$ and $j$ at a relative distance $d$ from each other, irrespective of their absolute spatial coordinates $y_i$ and $y_j$. Because the absolute coordinates $(y_i,y_j)$ of the pixel intensity pair $i$ and $j$ are no longer recorded, the GLCM is not a one-to-one mapping of the original gradient (unlike the 2D gradient map). For instance, in Fig. 3C, the pixel intensities $i=7$ and $j=0$ are located at a distance $d=(0,1)$, and they are represented in the 2D gradient map by a value of “1” at spatial coordinates $(y_i=7,y_j=8)$, as shown in Figs. 3A and 3B. However, the GLCM represents the same pair as an entry at $(i=7,j=0)$ as it remaps all 2D gradient map entries from Figs. 3C and 3G using a modulo $N_g$ operation. For example, the spatial coordinates $(y_i=7,y_j=0)$ from Fig. 3C are mapped modulo $N_g+1=9$ to GLCM coordinates $(y_i=7,y_j=0)$, which correspond to gray levels $(i=7,j=0)$ in GLCM. Although the $N_g\times N_g$ GLCM array can no longer be mapped back to the original image, it retains essential second-order spatial correlations of gray level intensities.

On the number of nonzero GLCM entries for a linear gradient

For any linear gradient $\mathrm{\nabla }$, the starting point of the GLCM has an index $i$ from the set $\{0,\mathrm{\nabla },2\mathrm{\nabla },\dots ,({\tilde{N}}_g-1)\mathrm{\nabla }\}$, where ${\tilde{N}}_g$ is the number of non-zero GLCM entries shown in Figs. 3D and 3H, i.e.:

(5) $${\tilde{N}}_g=1+\left[\frac{N_g-1}{|\mathrm{\nabla }|}\right].$$

In the above formula, $[\dots ]$ denotes the integer part. Each endpoint index $j$ of the GLCM is also expressed as $j=i+\mathrm{\nabla }$. This relationship indicates that the increment of endpoint indices in the GLCM, represented by $\mathrm{\Delta }j$, is equivalent to that of the starting point indices, denoted as $\mathrm{\Delta }i$, meaning that $\mathrm{\Delta }j=\mathrm{\Delta }i=\mathrm{\nabla }$, as illustrated by the horizontal and vertical double arrows in Fig. 3C for $\mathrm{\nabla }=1$ and in Fig. 3G for $\mathrm{\nabla }=2$. For example, in Fig. 3D and $\mathrm{\nabla }=1$ gray level per pixel along with Eq. (5) determines how many non-zero GLMC entries ${\tilde{N}}_g$ result from sampling the $N_g=8$ gray levels of the image, i.e., ${\tilde{N}}_g=1+[(8-1)/1]=8$. Similarly, for Fig. 3H with $\mathrm{\nabla }=2$ gray levels per pixel in conjunction with Eq. (5), it yields ${\tilde{N}}_g=1+[(8-1)/2]=4$.

One can notice that Fig. 3 displays the GLCMs for positive gradients $\mathrm{\nabla }>0$ and positive displacement vectors, such as $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)=(0,1)$. Reversing the direction of the gradient would merely shift all non-zero entries in the two-dimensional representation shown in Figs. 3C and 3G below the primary diagonal at a distance $i-j=d\cdot \mathrm{\nabla }<0$.

Marginal distribution of gray level differences $p_{x-y}$ for linear gradients

The marginal probability distribution $p_{x-y}(k)$, defined by Eq. (3) and visually represented in Fig. 2B, accounts for the sum of GLCM entries with specified gray level differences $k=j-i$. As observed in Figs. 3D and 3H, the lines parallel to the primary diagonal of the GLCM convey information about image gradients and represent the lines of constant gray level differences $p_{x-y}(k)$. For example, the GLCM primary diagonal entries have zero gray level differences, i.e., $k=j-i=0.$ Consequently, the sum of the primary diagonal elements, i.e., $p_{x-y}(0)={\sum }_{i=0}^{N_g-1}{\sum }_{j=0,|i-j|=0}^{N_g-1}p(i,j)$, is a zero gradient line because the gray level differences, i.e., the difference between the gray level value $i$ of the start (reference) point $(x_i,y_i)$ and the endpoint gray level intensity $j$ at $(y_j,x_j)$ along the displacement vector $d=(\mathrm{\Delta }x,\mathrm{\Delta }y)$, is $k=|i-j|=0$. Figure 3D illustrates that the GLCM of a gradient $\mathrm{\nabla }=1$ gray level per pixel along the vertical unit displacement vector $d=(0,1)$ contains all entries (except one) aligned along a parallel line with the primary diagonal at gray level differences $k=j-i=d\cdot \mathrm{\nabla }=1$. The sole exception is the GLCM entry at the discontinuity between the first period and the subsequent gradient repeats (see Fig. 3A and Fig. 3E). For example, the first period of the gradient in Fig. 3E and Fig. 3F ends with a gray level of $i=6$ in an image with $N_g=8$ gray levels and a gradient intensity $\mathrm{\nabla }=2.$ Therefore, its pair must have an intensity $j=i+\mathrm{\nabla }=8$, which is mapped modulo $N_g$ to $j=0$. It corresponds to $P(6,0)=1$ (remember that the wrapping around of GLCM in gray level spaces is done modulo $N_g$ because the gray level indices start at zero while the spatial coordinates wrap around with modulo $N_g+1$ operation because they begin at index 1). Since accounting for another period of the same gradient increases all nonzero entries of GLCM by one unit, from this point forward, one only calculates the GLCM for a single period of the gradient. To compute Haralick’s features, one uses the symmetry of the GLCM induced by periodic linear gradients such as those shown in Fig. 3.

One can observe from Fig. 3 that the nonzero GLCM entries parallel to the primary diagonal for a given gray level difference $k=j-i=d\cdot \mathrm{\nabla }$ begin at a distance of $d\cdot \mathrm{\nabla }$ from the first GLCM entry $p(0,0)$. The line of constant gray level differences $k=j-i=d\cdot \mathrm{\nabla }$ (dotted line parallel to the primary diagonal of the GLCM in Figs. 3D and 3H) starts at $p(0,d\cdot \mathrm{\nabla })$ and ends at $p(i=(m_1-1)\mathrm{\nabla },j=i+d\cdot \mathrm{\nabla })$, where $m_1$ is the number of GLCM entries along the gray level differences line with $k=j-i=d\cdot \mathrm{\nabla }$, which is:

(6) $$m_1={\tilde{N}}_g-|d|.$$

In the example depicted in Figs. 3A–3D, the GLCM for a unit vertical displacement $d=(0,1)$ in an image exhibiting a linear gradient of $\mathrm{\nabla }=1$ gray level per pixel and $N_g=8$ gray levels has a total of ${\tilde{N}}_g=8$ nonzero entries (from Eq. (5)), of which $m_1=7$ (see Eq. (6)) along the line of constant gray level differences $k=j-i=d\cdot \mathrm{\nabla }=1$. This line starts at $p(0,d\cdot \mathrm{\nabla })=p(0,1)$ and ends at $p(i=(m_1-1)\mathrm{\nabla },j=i+d\cdot \mathrm{\nabla })=p(6,7)$. Similarly, for the example shown in Figs. 3E–3H, $\mathrm{\nabla }=2$ gray levels per pixel and $N_g=8$, one gets a total number of GLMC entries of ${\tilde{N}}_g=4$ (from Eq. (5)), of which $m_1=3$ (see Eq. (6)) along the line of constant gray level differences $k=j-i=d\cdot \mathrm{\nabla }=2$ that starts at $p(0,d\cdot \mathrm{\nabla })=p(0,2)$ and end at $p(i=(m_1-1)\mathrm{\nabla },j=i+d\cdot \mathrm{\nabla })=p(4,6)$.

The GLCM always has exactly ${\tilde{N}}_g$ nonzero entries according to Eq. (5), of which, according to Eq. (6), $m_1$ are on the constant gray level differences line $k=j-i=d\cdot \mathrm{\nabla }$. The remaining $m_2$ nonzero GLCM entries have the endpoint coordinate $j$ always beginning at zero due to the wrapping around modulo $N_g$ in the gray level intensity space:

(7) $$m_2={\tilde{N}}_g-m_1=|d|.$$

Such GLCM entries are $p(i=m_1\mathrm{\nabla },j=0)$, $p(i=(m_1+1)\mathrm{\nabla },j=\mathrm{\nabla })$, and so on. One notices that all these new $m_2=|d|$ GLCM entries align along the line of constant gray level differences $k=j-i=-m_1\mathrm{\nabla }$, as shown in Figs. 3D and 3H. To summarize, the (unnormalized) marginal distribution of gray level differences $p_{x-y}$ for linear gradients represents the frequency of various combinations of pixel intensities that yield a specific absolute difference value $k=|j-i|$:

(8) $$p_{x-y}(k)=\{\begin{array}{cc}1,\hfill & \mathrm{f}\mathrm{o}\mathrm{r}k=j-i=d\cdot \mathrm{\nabla }\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}i=\{0,1,\dots ,m_1-1\}\mathrm{\nabla },\hfill \\ 1,\hfill & \mathrm{f}\mathrm{o}\mathrm{r}k=j-i=-m_1\mathrm{\nabla }\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}i=\{m_1,m_1+1,\dots {\tilde{N}}_g-m_1\}\mathrm{\nabla },\hfill \\ 0,\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}$$

Marginal distribution of gray level sums ${\mathit{p}}_{\mathit{x}\mathbf{+}\mathit{y}}$ for linear gradients

The previous section demonstrated that linear gradients are naturally represented by non-zero entries parallel to the primary diagonal of the GLCM. Thus, the marginal distribution of gray level difference $p_{x-y}$ arises naturally from GLCM symmetry. Other Haralick features require calculating the marginal distribution $p_{x+y}(s)$ for a given sum of gray level intensity $s=i+j$, where $s=\{0,1,\dots ,2(N_g-1)\}.$ One can utilize the GLCM symmetries caused by linear gradients and the corresponding marginal distribution $p_{x-y}(k)$ where $k=j-i=d\cdot \mathrm{\nabla }$ to streamline the calculation of the other marginal distribution $p_{x+y}$. Indeed, from $p_{x-y}(k)$, the $m_1$ nonzero endpoint gray level intensity are $j=i+k=i+d\cdot \mathrm{\nabla }$ where $i=\{0,\mathrm{\nabla },\dots ,(m_1-1)\mathrm{\nabla }\}$. Therefore, the elements of the marginal distribution $p_{x+y}(s)$ are $s=i+j=2i+d\cdot \mathrm{\nabla }$ with $i=\{0,\mathrm{\nabla },\dots ,(m_1-1)\mathrm{\nabla }\}$. Similarly, the second line of constant gray level differences is $k=j-i=-m_1\mathrm{\nabla }$ where $j=i-m_1\mathrm{\nabla }$ and $i=\{m_1\mathrm{\nabla },(m_1+1)\mathrm{\nabla },\dots \}$, which determines the marginal distribution $p_{x+y}(s)$ with $s=i+j=2i-m_1\mathrm{\nabla }$. In summary, the (un-normalized) marginal distribution of gray level sums $p_{x+y}$ for linear gradients indicates the frequency of various combinations of pixel intensities that total a specific value $s=j+i$:

(9) $$p_{x+y}(s)=\{\begin{array}{cc}1,\hfill & \mathrm{f}\mathrm{o}\mathrm{r}s=2i+d\cdot \mathrm{\nabla }\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}i=\{0,1,\dots ,m_1-1\}\mathrm{\nabla },\hfill \\ 1,\hfill & \mathrm{f}\mathrm{o}\mathrm{r}s=2i-m_1\mathrm{\nabla }\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}i=\{m_1,m_1+1,\dots {\tilde{N}}_{\mathrm{g}}-{\mathrm{m}}_1\}\mathrm{\nabla },\hfill \\ 0,\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}$$

Analytic scaling laws for Haralick features of linear gradients Comparison with numerical results

The previous subsection includes all the elements needed to estimate analytically any Haralick feature. In the following subsections, we derive analytical formulas for SA, SV, difference variance (DV), and Entropy based on the GLCMs symmetries derived in the previous subsections. Anticipating the results from the following subsections, the analytic scaling laws for Haralick features take the general form

$$f\propto N_g^{\alpha }|d{|}^{\beta }{\mathrm{\nabla }}^{\gamma },$$where the scaling exponents $\alpha ,\beta$ and $\gamma$ are derived from the GLCM symmetries as we will prove below.

To validate our theoretically predicted scaling laws for Haralick features, we performed numerical calculations using synthetic (computer-generated) gradient images. The predictions are represented by continuous lines in Figs. 4 and 5. At the same time, the corresponding numerical simulation results—based on the synthetic images described in “Synthetic Gradient Images”—are shown as discrete points with different symbols, as indicated in the figure legends.

To reduce plot clutter in Figs. 4 and 5, we present results only for odd intensity gradient values of $\mathrm{\nabla }\in \{1,3,5,7\}$ gray levels per pixels. In Fig. 4 the displacement vector magnitude was fixed at $|d|=1$, while the number of gray levels varied as $N_g\in \{16,32,64,128,256\}$. Conversely, in Fig. 5 the bit depth was set to $b=8$ bits ( $N_g=256$), while the vertical displacement vector magnitude varied as $|d|=1,\dots ,8.$ For each synthetic image with a given bit depth $b$ and gradient intensity $\mathrm{\nabla }$, we computed the GLCMs for each vertical displacement vector $d$ using Matlab’s $graycomatrix()$ function. For instance, the GLCM shown in Fig. 1C was obtained using $graycomatrix$(img,‘Offset’,[0 1], ‘NumLevels’, 4, ‘GrayLimits’, [], ‘Symmetric’,false). Additionally, when calculating all Haralick features, we consistently set the ‘Symmetric’ flag in graycomatrix() to true. Subsequently, we computed Haralick features from GLCMs using Matlab function $graycoprops()$.

For a single period of a linear gradient $\mathrm{\nabla }$ (see Fig. 3) all ${\tilde{N}}_g$ nonzero entries of the GLCM given by Eq. (5) have equal weight and are only aligned to two parallel lines to the primary diagonal as in Figs. 3D and 3H.

Sum average $f_6$

The SA indicates the uniformity of intensity values across the image texture. A higher SA value represents an even distribution of intensity sums between neighboring pixels. SA is defined as:

(10) $$f_6=\sum _{k=0}^{2(N_g-1)}k{p}_{x+y}(k).$$

A high SA implies that most pixel pairs have similar intensity sums, indicating a relatively uniform texture. A low SA suggests more significant variation in intensity sums between neighboring pixels, signifying a more textured appearance. From Eq. (10) with Eq. (9)

(11) $$\begin{array}{rl}& f_6=\frac{1}{{\tilde{N}}_g}(\stackrel{k=j-i=d\mathrm{\nabla }}{\overbrace{d\mathrm{\nabla }+(d\mathrm{\nabla }+2\mathrm{\nabla })+\dots (d\mathrm{\nabla }+2(m_1-1)\mathrm{\nabla })+}}\\ & \underset{k=j-i=-m_1\mathrm{\nabla }}{\underbrace{m_1\mathrm{\nabla }+(m_1\mathrm{\nabla }+2\mathrm{\nabla })+\dots (m_1\mathrm{\nabla }+2({\tilde{N}}_g-m_1)\mathrm{\nabla })}}=({\tilde{N}}_g-1)\mathrm{\nabla }=\mathrm{\nabla }\left[\frac{N_g-1}{\mathrm{\nabla }}\right].\end{array}$$

To simplify the calculation of $f_6$ above, we separated the contributions of the GLCM entries that are parallel to its primary diagonal at a distance of $k=j-i=d\mathrm{\nabla }$ from those on the line where $k=j-i=-m_1\mathrm{\nabla }$. Each of the two terms in Eq. (11) is an arithmetic series with the sum ${\sum }_{q=0}^{Q}a+2\mathrm{\nabla }q=a(Q+1)+Q(Q+1)\mathrm{\nabla }$. For $k=j-i=d\mathrm{\nabla }$ in Eq. (11) one uses $a=d\mathrm{\nabla }$ and $Q=m_1-1$ for $k=j-i=-m_1\mathrm{\nabla }$ one substitute $a=m_1\mathrm{\nabla }$ and $Q={\tilde{N}}_g-m_1$.

The first observation is that the theoretically predicted SA value given by Eq. (11) is independent of the gradient intensity $\mathrm{\nabla }$ (see the continuous lines in Fig. 4A) and the displacement vector $d$ (see the continuous lines in Fig. 5A) as summarized also in Table 1. Numerically computed Haralick feature SA confirms that its values are independent of gray level intensity gradients $\mathrm{\nabla }$ and increases linearly with the number of gray levels $N_g$ as shown in Fig. 4A. The exact formula in Eq. (11), which involves the discontinuous integer part function $[\dots ]$, is challenging to work with; however, by dropping the integer part operation, one finds a continuous approximate value ${\tilde{f}}_6\approx N_g-1$. This approximation demonstrates that $f_6$ scales linearly with $N_g$ (see the continuous lines Fig. 4A, which is also confirmed numerically by the linear increase of Haralick features with the number $N_g$ of gray levels shown in Fig. 4A. The second observation is that, numerical simulations shown in Fig. 5A confirm our theoretical prediction based on Eq. (11) that SA feature is independent of the displacement vector magnitude $|d|$. One notices, a slight error in approximating $f_6$ with ${\tilde{f}}_6$. For example, a $N_g=256$ gray level image and a gradient $\mathrm{\nabla }=2$ gray levels per pixel gives $f_6=2\left[(256-1)/2\right]=254$, which is slightly less than the simplified approximation ${\tilde{f}}_6=255$, but the error is under 0.4%. Even for gradients as large as $\mathrm{\nabla }=10$ gray levels per pixel, the error of approximating $f_6$ with ${\tilde{f}}_6\approx N_g-1$ is below 1%. This slight disagreement between the theoretical predicted SA value from Eq. (11) and the numerically computed values is emphasized in Fig. 5A. One can conclude that the gradient $\mathrm{\nabla }$ slightly decreases the SA value $f_6$, but the correction is negligible for small gradients $\mathrm{\nabla }<10$ gray levels per pixel. This fact is marked by the general attribute “independent” with an asterisk in Table 1.

Sum variance

The sum variance feature is defined as follows:

(12) $$f_7=\sum _{k=0}^{2(N_g-1)}{(k-f_6)}^2{p}_{x+y}(k),$$and can be analytically estimated for GLCM of linear gradients using the same strategy described above when deriving explicit analytical expression for SA in Eq. (11).

(13) $$\begin{array}{rl}& f_7=\frac{1}{{\tilde{N}}_g}(\stackrel{k=j-i=d\mathrm{\nabla }}{\overbrace{{(d\mathrm{\nabla }-f_6)}^2+{(d\mathrm{\nabla }+2\mathrm{\nabla }-f_6)}^2+\dots {(d\mathrm{\nabla }+2(m_1-1)\mathrm{\nabla }-f_6)}^2+}}\\ & \underset{k=j-i=-m_1\mathrm{\nabla }}{\underbrace{{(m_1\mathrm{\nabla }-f_6)}^2+{(m_1\mathrm{\nabla }+2\mathrm{\nabla }-f_6)}^2+\dots {(m_1\mathrm{\nabla }+2({\tilde{N}}_g-m_1)\mathrm{\nabla }-f_6)}^2}})=\\ & {\mathrm{\nabla }}^2({\tilde{N}}_g^2/3-{\tilde{N}}_{g}d+d^2-1/3).\end{array}$$

To accurately predict the scaling law of SV features from the exact formula given by Eq. (13), one could eliminate the integer part function from the definition of ${\tilde{N}}_g$ and utilize an approximate estimate:

(14) $${\tilde{f}}_7=((N_g-1{)}^2/3-(N_g-1)d\mathrm{\nabla }+d^2{\mathrm{\nabla }}^2-1/3{\mathrm{\nabla }}^2).$$

The discrepancy between the true $f_7$ (Eq. (13)) and the approximate estimate ${\tilde{f}}_7$ is minor but can reach several percentage points. For example, the largerst error occurs for $N_g=256$, $\mathrm{\nabla }=7$, and $|d|=8$, which is approximately 3.33%.

Based on Eq. (14), one notice that SV scales quadratically with $N_g$. Indeed, the second term in Eq. (14), which is linear in $N_g$, is always smaller than the first term, quadratic in $N_g$ if $d\mathrm{\nabla }<N_g$. This condition is fulfilled because the product $d$ pixels times $\mathrm{\nabla }$ gray levels per pixel is the number of gray levels variation across an image, which cannot be larger than $N_g$. Numerical simulations confirmed our analytical prediction of a quadratic scaling law for $f_7$ with $N_g$, as shown in Fig. 4B. One also notices from Fig. 4B that for fixed displacement vector magnitude $|d|$, numerical values of SV are independent of gradient intensity $\mathrm{\nabla }$ as predicted analytically by Eq. (14).

For an image with a fixed number of gray levels $N_g$ and gradient $\mathrm{\nabla }$ gray levels per pixel, the second term in Eq. (14) dominated SV’s dependence on $|d|$. This is because $(N_g-1)d\mathrm{\nabla }>d^2{\mathrm{\nabla }}^2$, which reduces to $N_g-1>d\mathrm{\nabla }$. This was shown above to be true for all images. Furthermore, the second term in Eq. (14) $(N_g-1)d\mathrm{\nabla }$ is also larger than the fourth term $1/3{\mathrm{\nabla }}^2$ because $(N_g-1)d>1/3\mathrm{\nabla }$ even for the smallest possible displacement vector with $|d|=1$. As a result, the linear term $|d|$ is the primary influence in the scaling law of $f_7$, which aligns with our numerical simulations shown in Fig. 5B. As noticed from Fig. 5B, for a fixed displacement vector magnitude $|d|$, SV linearly changes with the gradient intensity $\mathrm{\nabla }$ as predicted analytically (see also Table 1).

Difference variance

The definition of difference variance is:

(15) $$f_{10}=\sum _{k=0}^{2(N_g-1)}{(k-DA)}^2{p}_{x-y}(k),$$where the DA is given by $DA={\sum }_{k=0}^{2(N_g-1)}k{p}_{x-y}(k)$. The evaluation of DA is straightforward and follows from Eq. (8) since all GLCM entries are equal weight:

$$DA=\frac{1}{{\tilde{N}}_g}(m_{1}d\mathrm{\nabla }+m_2(-m_1\mathrm{\nabla }))=0.$$

As a result, the DA reduces to

(16) $$f_{10}=\sum _{k=0}^{2(N_g-1)}{k}^2{p}_{x-y}(k)=\frac{1}{{\tilde{N}}_g}(m_1(d\mathrm{\nabla }{)}^2+m_2(-m_1\mathrm{\nabla }{)}^2)={\mathrm{\nabla }}^{2}d({\tilde{N}}_g+1-d).$$

To infer the asymptotic scaling law exponents from the exact formula of DS given by Eq. (16), one drops the integer part function from ${\tilde{N}}_g$ and uses an approximate formula ${\tilde{f}}_{10}\approx \mathrm{\nabla }d(N_g-1)+(1-d){\mathrm{\nabla }}^2\approx \mathrm{\nabla }d{N}_g$, which suggests the scaling law

$$f_{10}\propto N_g|d|\mathrm{\nabla }.$$

The theoretically predicted linear scaling with $N_g$ is confirmed by numerical simulations shown in Fig. 4C, for a fixed $|d|=1$ pixel and slopes that increase linearly with the gradient intensity $\mathrm{\nabla }$.

The scaling of experimental $f_{10}$ with the displacement vector $d$ exhibits a linear dependence with a slope proportional to the gradient $\mathrm{\nabla }$. Additionally, the plot of the theoretical prediction from Eq. (16) shows some deviation from linearity for large gradients. This is expected because ${\tilde{f}}_{10}$ neglects the contribution of the term $(1-d){\mathrm{\nabla }}^2$ compared to $\mathrm{\nabla }d(N_g-1)$. However, the contribution of the neglected term increases quadratically with the gradient $\mathrm{\nabla }$ and could become significant for images with large gradients (see Table 1).

Entropy

The Haralick features discussed thus far are derived from different moments of the marginal distribution of either the difference intensity (see Eq. (8)) or the sum intensity (see Eq. (9)). In contrast, entropy employs a logarithmic scale to compute features from the GLCM. The definition of the entropy feature is:

(17) $$f_9=-\sum _{i=0}^{N_g-1}\sum _{j=0}^{N_g-1}p(i,j)\log (p(i,j)).$$

Entropy reaches its maximum value when a probability distribution is uniform (entirely random texture) and its minimum value of 0 when all grayscale values in the image are the same. If the entropy $f_9$ is defined using the base-2 logarithm ${\log }_2()$, then $f_9$ is measured in bits. While we only examined the entropy feature $f_9$, employing the marginal distributions of pixel intensity sums from Eq. (9) along with the detailed calculation examples for the SA and SV features, one can easily deduce the scaling law of SE defined by:

(18) $$f_8=-\sum _{k=0}^{2(N_g-1)}{p}_{x+y}(k)\log (p_{x+y}(k)).$$

Similarly, by using the marginal distribution of pixel intensity difference from Eq. (8) and the detailed calculation examples provided for the DA feature, one can easily derive the scaling law of DE defined by:

(19) $$f_{11}=-\sum _{k=0}^{N_g-1}{p}_{x-y}(k)\log (p_{x-y}(k)).$$

Calculating entropy is straightforward because all GLCM entries carry equal weight, leading to:

(20) $$f_9=-\sum _{i=0}^{N_g-1}\sum _{j=0}^{N_g-1}p(i,j)\log (p(i,j))=\log ({\tilde{N}}_g).$$

As seen from the numerical simulation results presented in Fig. 4D, the theoretical scaling law derived from Eq. (20) captures the general logarithmic trend of the entropy. However, it slightly underestimates it (see Table 1). Numerical simulations illustrated in Fig. 5D confirm that Entropy feature $f_9$ is independent of the magnitude of the displacement vector, as predicted by Eq. (20), and also slightly underestimates the actual values. The discrepancy arises from an offset constant $\epsilon$ used in estimating the entropy from images where $\log (p(i,j)+\epsilon )$ was employed instead of $\log (p(i,j))$ to prevent the entropy singularity for sparse GLCM.