1 Introduction

Texture is a fundamental property of images that describes the spatial arrangement of intensities or colors. It plays a crucial role in various scientific fields, including object recognition, medical image analysis, and product quality control. Texture classification is essential for tasks such as identifying objects in images, detecting abnormalities in medical scans, and ensuring the consistency of manufactured products. A previous step for this purpose is the accurate characterization of image features, which is precisely the aim of this paper. Traditional methods for extracting and analyzing textures can be computationally expensive and susceptible to image quality issues such as noise, blur, and color variations.

To overcome these limitations, there is a growing need for robust and efficient methods able to accurately discriminate textures from various sources, including natural product images such as those of rocks, fabrics, and skin, which exhibit complex and diverse textures. A successful texture discrimination method should be able to effectively distinguish between them, even under varying image quality conditions. Moreover, rotational invariance is a desirable feature in algorithms. This means that the algorithm can accurately identify images regardless of their orientation or rotation [19]. Thus, a reliable method should account both for image rotation and other types of rigid transformations.

The development of such methods has the potential to revolutionize various scientific domains. In object recognition, texture identification can aid in discriminating objects based on their unique textural signatures. In medical image analysis, it can help in detecting abnormalities and assessing tissue health by analyzing the textures of diseased tissues. In product quality control, it can be used to ensure consistent product quality by identifying deviations in texture patterns. In the fashion industry, it could be a valuable tool for identifying and understanding consumer preferences for different types of garments.

The use of information-theory-related methods in one dimension is not new in the analysis of time series. In fact, it has been successfully used for the analysis of dynamical systems in different scientific domains. For example, [36] uses wavelet entropy to characterize electroencephalography records across different stages of brain activities; [28] proposes the use of permutation entropy and statistical complexity for selecting the most convenient randomizing technique; in a similar vein, [38] explores the relationship between entropy and epochs of rapid climate change in the Holoce, [37] highlights the robustness of permutation entropy and statistical complexity quantifiers to characterize a long list of nonlinear chaotic maps. [42] finds that permutation entropy and complexity reveal key dynamical features of a semiconductor laser, and [5] unveils interest rate manipulation using global and local information theory quantifiers.

There are several studies that use entropic measures for studying images. For example, [45] employs permutation entropy and statistical complexity to examine the evolution of internet memes, particularly focusing on changes in visual complexity over time. Zhang et al. [47] adopt a bottom-up approach, developing predictive models based on features extracted from paintings. Their work reveal a correlation between key painting features and permutation entropy and statistical complexity. In [18] permutation entropy (PE) and amplitude-aware permutation entropy (AAPE) in two dimensions are used to analyze synthetic and biomedical images. They find that both PE and AAPE were effective in differentiating textures and identifying irregularities. In a subsequent study, [17] extract quantitative information from synthetic and biomedical images using ensemble dispersion entropy, ensemble permutation entropy, and ensemble sample entropy. Their results demonstrate the ability of ensemble approaches to detect randomness, periodicity, and tissue types in the analyzed images. Regarding image assessment, [20] highlights the significance of entropy as a measure of randomness and information content in capturing aesthetic judgments. The authors argue that entropy reflects the complexity and unpredictability of an image, contributing to its overall aesthetic evaluation. More recently, [32] proposes the Robust Directional Median Pattern to characterize textures, resulting in more competitive than other local binary pattern methods for texture classification.

However, as [13] recalls, there is a need to study spatial structure and pattern in more than one dimension. In particular, it is not trivial to understand how information is assembled across a two-dimensional lattice, and how to extract such information. By ‘global’ and ‘local’ quantifiers we mainly refer to the Shannon entropy and the Fisher information measure, respectively. Motivated by this fact, Ribeiro and coauthors [33], introduce an extension of the complexity-entropy causality plane to evaluate the complexity of two-dimensional patterns, as are fractal images, liquid crystal textures and Ising surfaces. Roughly speaking, they consider an image as a \(m\times n\) matrix that is split into \(D_m\! \times \! D_n\) submatrices. Each submatrix is converted, following a row-wise order, into a symbol and then, they derive a time series from the calculated symbols and compute the permutation information theory quantifiers. Finally, they use the complexity-entropy causality plane to evaluate the complexity of image patterns. Later, Zunino and Ribeiro [48], apply this extension to study texture discrimination in real images. Sigaki et al. [41] also apply this extended method to analyze different styles in historic art paintings and classify their artworks. It is important to mention that these works approach the topic of higher dimensions only for the complexity quantifier, which is interpreted, as said before, as a global one. Up to our knowledge, there is no proposed method extended to higher dimension arrangements to calculate the Fisher Information Measure, neither the corresponding Fisher-entropy plane. The value of the Fisher Information Measure defined in Eq. 8, for a discrete probability distribution \(P\!=\!\big \{p_j, j\!\in \!\{1,\ldots ,N\}\big \}\), strongly depends on the j-ordering –hence its condition as a local quantifier. The key point is that there is no a uniquely predetermined order, when P is defined on two –or higher– dimension patterns, to convert a matrix into a symbol.

In this sense, this paper presents a novel two-step method that first transforms a two-dimensional image into a one-dimensional time series using the Hilbert curve and then computes three information theory quantifiers using the Bandt & Pompe symbolization method. The method extracts the intrinsic characteristics of the images and allows their characterization and discrimination. The election of the planar Hilbert curve arises from the need to convert, as best as possible, spatial correlations into temporal ones.

We first evaluate the proposed method on images generated by self-similar multifractal surfaces and by Brownian surfaces, demonstrating, in the last case, its ability to accurately separate textures based on their Hurst exponents. Next, we apply the method to the Brodatz image database, which has been widely used in texture studies, and demonstrate its robustness to image rotation, mirror-symmetry, and color variations.

Our contribution is based on the development of a novel method for analyzing images by transforming them into a regular time series that offers several advantages over existing approaches: (i) It eliminates directional biases inherent in row/column-wise scan, ensuring a more balanced representation of the image’s texture; (ii) It leverages both the Hilbert curve’s space-filling and non-privileged direction properties to generate a path that visits each pixel once, ensuring no data redundancy; (iii) it utilizes information theory quantifiers, permutation entropy, permutation statistical complexity and Fisher information measure, to provide a robust and informative characterization of different textures; (iv) it results robust under rigid images transformations; (v) it is also invariant to image scaling; (vi) it can be generalized in a natural way to more than two dimensions –as the Hilbert curve is defined in \({\mathbb {R}}^n\)– and be effectively applied to discriminate n-dimensional patterns.

The remainder of the paper is organized as follows. Section 2 provides a summary of the information quantifiers considered, followed by Sect. 3 where the proposed method is presented. Section 4 discusses the key findings and finally Sect. 5 summarizes the paper’s contributions.

2 Statistical information quantifiers

This section describes the information theory quantifiers that will be used in the empirical section of the paper.

2.1 Shannon entropy

In the context of discrete random variables, we consider an arbitrary variable denoted as X. Our primary objective is to quantify the information content obtained upon observing this variable. The quantification of this information relies upon the associated probability distribution, denoted as p(x). In the field of dynamical systems, Shannon entropy emerges as a metric capturing the smoothness of the probability distribution p(x). It inherently characterizes the degree of disorder within the system, with respect to the likelihood of various accessible states. In situations where an event is certain (\(p(x)\!=\!1\)), the information gained is minimal, as the certainty of its occurrence was already established. Conversely, in the case of improbable events (\(p(x)\!=\!0\)), the revelation of such occurrences imparts a substantial amount of information. Indeed, the highest level of uncertainty regarding an event manifests itself when the probability distribution is uniform. Following the work of Claude Shannon [40], we can define Shannon’s Entropy as:

$$\begin{aligned} S[P]=-\sum _{j=1}^N p_j \ln (p_j) \end{aligned}$$
(1)

Where \(P\!=\!\big \{p_j, j\!\in \! \{1,\dots , N\}\big \}\) represents a specific probability distribution.

Shannon’s Entropy serves as a ‘global’ measure and exhibits limited sensitivity to localized perturbations within the distribution. To mitigate this limitation and facilitate comparative analysis, Shannon’s Entropy can be normalized to a range of values between 0 and 1. This normalization is achieved by dividing the entropy by its maximum attainable value \(S_{\max }\), and defined as:

$$\begin{aligned} \mathcal{H}[P]=\frac{S[P]}{S_{\max }}=\frac{-\sum _{j=1}^N p_j \ln (p_j)}{\ln N} \end{aligned}$$
(2)

2.2 Statistical complexity

Analyzing time series using only Shannon entropy may not provide a complete understanding of a system’s underlying dynamics. As noted by [12], entropy alone is insufficient to capture the intricate organization and structure that often define complex systems. To gain a more comprehensive perspective, statistical complexity measures are introduced. These measures can reveal hidden patterns and correlations within time series, enabling a more accurate system behavior characterization. Time series with a rich interplay of patterns are assigned higher complexity values, reflecting a greater degree of structure and organization. In contrast, highly ordered or random sequences, such as perfectly periodic or entirely unpredictable patterns, are assigned lower complexity values. This ability to differentiate systems based on their level of internal organization underscores the utility of statistical complexity measures.

Based on the functional form defined by [22, 24] introduces a complexity measure that reads:

$$\begin{aligned} \mathcal{C}_{JS}[P]=\mathcal{Q}_J[P,P_e] \,\mathcal{H}[P] \end{aligned}$$
(3)

where \(P_e\) is the uniform distribution (i.e., the equiprobability given by \(p_j\!=\!1/N\) \(\forall j\!\in \!\{1,\dots , N\}\)), and \(\mathcal{Q}_J\) is the disequilibrium between the observed and the \(P_e\) distribution, defined in terms of the Jensen-Shannon divergence, which could be regarded as statistical distance between two points of an N-dimensional probability space [23, 27]. Let consider:

$$\begin{aligned} \mathcal{Q}_{J}[P,P_e]=\mathcal{Q}_0 \,\mathcal{J}[P,P_e] \end{aligned}$$
(4)

with:

$$\begin{aligned} \mathcal{J}[P,P_e]=S[(P\!+\!P_e)/2]-S[P]/2-S[P_e]/2 \end{aligned}$$
(5)

where \(\mathcal{Q}_0\) is a normalization constant that represents perfect certainty (i.e. one of the probability states is equal to one and the remaining \(p_j\) are equal to zero), namely:

$$\begin{aligned} \mathcal{Q}_0=-2\left[ \frac{N\!+\!1}{N}\, \ln (N\!+\!1) - 2 \ln (2N)+\ln (N) \right] ^{-1} \end{aligned}$$
(6)

The statistical complexity outlined above is calculated by multiplying two normalized measures of uncertainty: Shannon entropy and Jensen-Shannon divergence. However, this complexity is not simply a function of one distribution’s uncertainty. Instead, it is influenced by the system’s deviation from a uniform reference state (\(P_e\)), highlighting the interplay between order and disorder in the system.

2.3 Fisher information measure

In a landmark contribution to statistical science, [15] introduces the concept that nowadays is known as Fisher Information Measure (FIM), a versatile metric that quantifies the amount of information contained within a set of data about an unknown parameter. Beyond its role in parameter estimation, FIM serves as a measure of the informational yield of a data set, reflecting the ‘quality’ of the measurements and providing insights into the underlying system or phenomenon [16].

Unlike Shannon entropy, FIM distinguishes itself by incorporating derivative terms, effectively gauging the sensitivity of the probability density function to reordering over the independent variable. This heightened sensitivity proves particularly valuable in scenarios where a clear ordering exists, such as in dynamic systems where time naturally serves as an ordering variable. Thus, FIM quantifies the ‘gradient content’ of a distribution and reads:

$$\begin{aligned} \mathcal{F}[f]=\int _{\Delta } \frac{1}{f(x)}\left[ \frac{df(x)}{dx}\right] ^2 dx = 4 \int _{\Delta } \left[ \frac{d \psi (x)}{dx} \right] ^2 \end{aligned}$$
(7)

The standard FIM definition involves dividing by the distribution itself, which becomes problematic when the distribution approaches zero. Using real probability amplitudes (\(\psi (x)\)) avoids this issue and shows that FIM simply measures the gradient in \(\psi (x)\) [15, 16]. This sensitivity to local variations makes FIM a ‘local’ quantifier.

For discrete systems with N states, discretizing FIM introduces complications like losing shift-invariance. We focus on a normalized discrete FIM \((0 \le {{\mathcal {F}}}\le 1)\) given by:

$$\begin{aligned} {{\mathcal {F}}}[P]~=~F_0~\sum _{i=1}^{N-1}\big [(p_{j+1})^{1/2} - (p_{j})^{1/2}\big ]^2 \end{aligned}$$
(8)

This discretization is considered well-behaved in discrete settings [43]. The normalization constant \(F_0\) reads:

$$\begin{aligned} F_0~=~\left\{ \begin{array}{cl} 1 & \quad \text{ if } \,p_{j^*}\!=\! 1\text{, } \text{ for } j^*\!=\! 1 \text{ or } j^*\!=\! N \text{ and } \,p_{j}\!=\!0 \ \ \forall \, j\!\ne \! j^* \\ 1/2 & \quad \text{ otherwise } \end{array} \right. \end{aligned}$$
(9)

A highly ordered system (most \(p_j\) near zero, one \(p_k\) near 1) has low Shannon entropy (\({{\mathcal {H}}}\!\sim \!0\)) and high normalized FIM (\({{\mathcal {F}}}\!\sim \!1\)). Conversely, a disordered system (all \(p_j\) similar) has high entropy (\({{\mathcal {H}}}\!\sim \!1\)) and low FIM (\({{\mathcal {F}}}\!\sim \!0\)). This discrete FIM behaves opposite to Shannon entropy, except for periodic motions [30, 31].

The Fisher Information Measure can be interpreted as capturing the roughness of a probability density function, quantified by differences between adjacent probability states. This requires establishing a specific order for these states. The relevant literature often adopts either the ordering proposed by [21] or the ‘Lehmer code’ described by [39], the latter being generated through an algorithm based on the factoradic number system. In this paper, we follow the Lehmer code ordering, as supported by previous findings in the literature [5, 30, 31].

2.4 Bandt and Pompe symbolization method

To calculate the aforementioned quantifiers, it is essential to establish the probability density distribution of the system under investigation. Several alternative methods are available for this purpose, including amplitude-based approaches [28] and binary symbolic dynamics [29], among others. However, two decades ago, Bandt and Pompe (BP) [4] introduce a powerful technique that has since been widely applied across various fields, including physics, biology, finance, and climate science. This method, which relies on the ordinal relations between neighboring values, converts the observed time series into symbolic sequences, enabling the extraction of meaningful features and patterns from complex temporal data.

The BP technique operates as follows. Let consider a time series \(\mathcal{S}(t)\!=\!\{x_t\}_{1\le t \le N}\), an embedding dimension \(D\!>\!1\), and an embedding delay \(\tau\), with \(D,\tau \!\in \!{{\mathbb {N}}}\). Then, we are in a position to define patterns of order D, generated by:

(10)

Then, patterns are evaluated according to the ordinal position of their components. For example, let \(\{1,6,3,2,8,9,5\}\) be a time series, \(D\!=\!4\) and \(\tau \!=\!1\) the embedding and delay dimensions. The following patterns are constructed: (1, 6, 3, 2), (6, 3, 2, 8), (3, 2, 8, 9), (2, 8, 9, 5), which are transformed into the following symbols by substituting the original values by the respective ordered position: (0, 3, 2, 1), (2, 1, 0, 3), (1, 0, 2, 3), (0, 3, 1, 2). Therefore, BP transforms the original time series into a series of symbols of length D. These symbols, by the nature of their construction, include temporal information, which grows along with D.

Given an embedding dimension D, there are D! patterns. The frequency of appearance of these patterns in the time series could be used as a straightforward measure to construct the probability density function associated to the symbolic time series \(P\!=\!\big \{p(\pi _i),\,i\!=\!1,\dots , D!\big \}\), where:

$$\begin{aligned} p(\pi _i)=\frac{\#\big \{s\!:\,s\le N\!-\!(D\!-\!1)\,\tau ;\ (s) \text{ has } \text{ type } \pi _i\big \}}{N\!-\!(D\!-\!1)\,\tau } \end{aligned}$$
(11)

It was shown in previous research that Brownian motion, fractional Brownian motion, and fractional Gaussian noise processes include all permutation patterns in their respective PDF. In other words, there are no ‘forbidden patterns’, and the relative frequency of each pattern will depend on the type of stochastic process under analysis. On the contrary, chaotic processes could include them, as some orbits are never observed due to their deterministic nature [2, 34].

The BP approach offers several advantages over conventional methods based on range partitioning. It is invariant to nonlinear monotonic transformations, meaning that nonlinear drifts or scalings introduced by measurement devices will not affect the estimation of quantifiers-an important property when dealing with experimental data. Additionally, the BP method is simple and computationally efficient, requiring only two parameters: the pattern length or embedding dimension (D) and the embedding delay (\(\tau\)). Furthermore, it can be applied to any time series, whether regular, chaotic, noisy, or derived from real-world data. The only condition for its applicability is a weak stationarity assumption, which requires the probability of \(x_t\!<\! x_{t+k}\) be independent of t for \(k\!=\!D\).

Overall, the BP method of symbolic analysis of time series has proven to be a versatile and robust tool for characterizing complex signals. It offers a unique perspective on the underlying dynamical system, gaining deeper insights into a wide range of phenomena.

2.5 Complexity entropy and entropy fisher causality planes

In a breakthrough paper, Rosso and coworkers [35] show that the joint study of the permutation entropy and permutation statistical complexity provides useful information to discriminate stochastic and chaotic processes. Causality indicates the temporal correlations between successive samples in the ordering suggested by Bandt and Pompe.

For each probability distribution, \(P\!=\!\{p(\pi )\}\), the pair \(\big (\mathcal{H}[P],\mathcal{C}_{JS}[P]\big )\) is represented in the plane \(\mathcal{H}\!\times \! \mathcal{C}\), so called Complexity Entropy Causality Plane (CECP). Although such pairs must belong to \([0,1]\!\times \![0,1]\), due to their normalized components, Martin, Plastino and Rosso [27] explain that in the CECP, all the possible values are bounded by two continuous curves representing the maximum and minimum values of \(\mathcal{C}_{JS}\) as a function of \(\mathcal H\). The CECP has been effectively used to rank stock markets according to their informational efficiency [49], to classify commodities [14], to unveil the chaotic dynamics of a semiconductor laser [42].

An additional insight into the hidden characteristics of dynamical systems could be provided by the so called Fisher Entropy Causality Plane (FECP), introduced by Vignat and Bercher [46], in which the pair \(\big (\mathcal{H}[P],\mathcal{F}[P]\big )\) is represented in the plane \(\mathcal{H}\!\times \!\mathcal{F}\). It was used, for example, to detect the change in the observed stochastic process due to manipulation of data [5].

Given the importance of both global and local perspectives, this article will classify textures using both planar representations.

3 Proposed method

It seems reasonable to assume that the CECP can be used to analyze properties of images, considering the value and relative arrangement of their pixels. The key problem is how to extend the ordering proposed by BP from 1D to 2D, or conversely, how to transform a 2D image into a 1D time series. Transforming a digitized surface into a 2D image necessitates a strategic reading approach to extract meaningful information. Conventional methods, such as reading by rows or simultaneously processing multiple rows, overlook the inherent spatial relationships between pixels, introducing potential biases due to preferential reading directions.

To address these limitations, we propose a novel methodology that traverses the entire image in a non-repeating, direction-diverse manner, ensuring that each pixel is visited once and only once. We achieve this by employing the Hilbert curve, a space-filling curve that generates a path that visits every pixel while maintaining a high degree of spatial locality.

3.1 The Hilbert curve

This well-known curve, is a fractal continuous curve in \({\mathbb {R}}^n\), with full Hausdorff dimension n, due to its space-filling property. The curve is obtained as a limit of polygonal curves of increasing lengths, called levels, as it is shown in Fig. 1 for \(n\!=\!2\). Notice that the higher the level the more space is filled by the polygonal curve. If the unit square is divided into \(2^N\!\!\times \!2^N\) sub-squares, then a Hilbert curve of level N is needed to scan all sub-squares.

A very important feature of space-filling curves is that they can map multidimensional data to one dimension while preserving locality of the data points [3]. This means that two close points in a one-dimensional arrangement correspond to two close points in the n-dimensional arrangement before unfolding (the converse is not true). This can be seen in Fig. 1 by noting that consecutive nodes in a polygonal curve, correspond to neighboring sub-squares.

Fig. 1
figure 1

Levels of planar Hilbert curve. a level 1, b level 2, c level 3, d level 6

This property is very convenient when the local correlation of the data is to be preserved. Additionally, while other curves can also unfold a two-dimensional arrangement with this property, the Hilbert curve does so more effectively. This efficiency comes from its self-similarity, which follows a recursive pattern, and it accomplishes this fact without favoring any particular direction of travel. These features make the Hilbert curve to be widely used in data-based applications, as multi-dimensional indexing methods, data structures, parallel computing, image processing, and encryption algorithm [8, 44].

The resulting time series, extracted from the Hilbert curve-based scan, captures the spatial characteristics of the image while mitigating the inherent biases of conventional methods. We then employ information theory quantifiers described in Sect. 2, using the CECP and FECP, to relate the time series features to distinct textures. These planar representations provide a comprehensive and insightful understanding of the texture characteristics.

4 Results

The proposed method is tested on simulated images and in a widely used texture database, namely: (i) self-similar multifractal surfaces, (ii) Brownian surfaces, (iii) black and white Brodatz dataset, and (iv) colored Brodatz dataset. The experiments not only include varying the value of the parameters in the calculation of the information quantifies, but also include the rotation and mirror-symmetry transformation of the images.

4.1 Self-similar multifractal measures as complex images

The complexity of an image is determined by the spatial distribution of a scalar quantity, as it could be the color or gray-level at every pixel. This quantity, when normalized, can be considered as a probability mass distribution, \(\mu\), supported on the image. The way \(\mu\) is spread on the region characterizes the type of image in question. When the mass density is so highly irregular that its concentration obeys different local power laws, and at the same time, the sets of points having the same local mass concentration define different fractals subsets on the support, then \(\mu\) is called a multifractal measure. It posses a very rich structure at every scale which yields to a whole set of fractal dimensions, usually studied by means of various classes of multifractal spectra [11].

In this section we investigate a multifractal mass distribution defined on the unit square as a planar image, by means of information-theory methods rather than by terms of a classical multifractal dimensional analysis.

4.1.1 Multinomial multiplicative cascades

An example of a process generating a self-similar multifractal measure (SMM), is a deterministic one, where an Euclidean set, \(E\!\subset \!{\mathbb {R}}^n\), is fragmented into smaller equal-size pieces according to a recursive fixed rule, and simultaneously, the total mass (\(\mu (E)\!=\!1\)) is distributed on the pieces following another non trivial recursive rule. This process is called a multiplicative cascade, and generates an SMM called a multinomial measure, which has features typical of a large class of multifractal ones [10].

To illustrate the process, let \(E\!=\![0,1]\!\subset \!{{\mathbb {R}}}\) and \(p_1,p_2\!>\!0\), with \(p_1+p_2\!=\!1\). A binomial SMM, \(\mu \!=\![p_1,p_2]\), is constructed on E, as follows: at first step, E is divided into two sub-intervals of length 1/2, and an initial uniform density is redistributed by assign mass \(p_1\) to the left sub-interval, and mass \(p_2\) to the right one. In the following steps, the equal-size subdivision of intervals is applied successively, the same as the mass on each sub-interval is rescaled by a factor of \(p_j\), as it is, left \(j\!=\!1\), or, right \(j\!=\!2\). The process is illustrated in Fig. 2. The last panel exhibits the extreme irregularity of measure \(\mu\).

Fig. 2
figure 2

Construction of a SMM \(\mu \!=\![p_1,p_2]\) (binomial cascade) on the unit interval with \(p_1+p_2\!=\!1\), \(p_j\!>\!0,j\!=\!1,2\), starting from the uniform distribution. A scheme of the first three steps and an advanced step of the recursive process are shown (see the text for the explanation)

4.1.2 Worked example

The analogous of the binomial SMM constructed on [0, 1], for the two-dimensional case, is a similar construction but on the unit square, \(E\!=\![0,1]\!\times \![0,1]\!\subset \!{{\mathbb {R}}}^2\). We will consider a quadrinomial SMM, \(\mu \!=\![p_1,p_2,p_3,p_4]\), obtained from the same recursive process described above, by dividing each square into four equal size sub-squares, with the values: \(p_1\!=\!0.2434,p_2\!=\!0.2522,p_3\!=\!0.2566\) and \(p_4\!=\!0.2478\) (\(\sum _{j=1}^4p_j\!=\!1\)) (see Fig. 3). The image obtained at the \(10^\textrm{th}\) step of the process and shown in Fig. 4(a), will be used as an example of a multifractal-type image. It consists of a \(2^{10}\!\times \!2^{10}\) arrange of \(\mu\)-values which will be scanned by a Hilbert curve of level 10 and then unfolded to get the corresponding time series.

Fig. 3
figure 3

Three steps of SMM construction, \(\mu \!=\![p_1,p_2,p_3,p_4]\) on the unit square, \(p_1\!=\!0.2434,p_2\!=\!0.2522,p_3\!=\!0.2566\) and \(p_4\!=\!0.2478\). a Measure values attached to sub-squares at step 1, b The same but at step 2, c The same but at step 3 (only a few). The colors are in accordance to the measure values

Fig. 4
figure 4

Multifractal-type image. a Measure values attached to sub-squares at step 10 (only the colors), b 3-dimensional rendering of (a) showing the SMM \(\mu\). The colors are in accordance to the measure values. (\(p_1\!<\!p_4\!<\!p_2\!<\!p_3\))

In order to analyze the performance of the proposed method on different types of images, we start by producing two different new ones, built from the same SMM surface of Fig. 4. The first one is obtained by reordering the \(\mu\)-values of the \(2^{10}\!\times \!2^{10}\)-arrange, from lowest to highest, yielding a complete regular image as shown in Fig. 5. The second one, by randomly reordering the same \(\mu\)-values, yields a fully noisy image as shown in Fig. 6. These last two images will be also scanned by a Hilbert curve of level 10 and subsequently unfolded into time series.

Fig. 5
figure 5

a Regular image obtained by ordering the \(\mu\)-values of the multifractal surface in Fig. 4 from lowest to highest, b 3-dimensional rendering of (a)

Fig. 6
figure 6

a Random image obtained by mixing randomly the \(\mu\)-values of the multifractal surface in Fig. 4, b 3-dimensional rendering of (a)

For these three time series given by the images, henceforth the multifractal, the ordered and the randomized cases, we compute the values of the corresponding information theory quantifiers according to the methods described in Sec. 2, for fixed embedding parameter \(D\!=\!6\) and delay parameter \(\tau \!=\!1\), to maintain the spatial resolution of the image as much as possible. We also consider the following rigid transformations on each of the original images: \(\frac{\pi }{2},\pi\) and \(\frac{3\pi }{2}\) angle rotations, and a mirror-symmetrical image, so a set of five images (and corresponding time series) are related to each case.

Fig. 7
figure 7

CECP (a) and FECP (b) for the SMM (red), the ordered (blue) and the randomized (green) cases. A close up on each group is made to show the true planar locations of the five points corresponding to the original, the rotated (\(\frac{\pi }{2},\pi\) and \(\frac{3\pi }{2}\)) and mirror-symmetrical images (color figure online)

Figures 7(a) and (b) display the CECP and FECP, respectively, for the fifteen data points. It can be appreciated the net separation in the planar location of the data corresponding to the multifractal, the ordered and the randomized cases, that shows the ability of the purposed method to discriminate the images type between the three classes. The invariance of the results under rotations and mirror-symmetry transformations is also fully evident, since the points of each group are located almost in the same place. It is necessary to make a close-up of each group to be able to distinguish between the five data points, and in the randomized case there are still points in the same place. For the multifractals, there is a very little discrepancy that separates two data points corresponding to the original and the \(\pi\) angle rotated images from three data points corresponding to \(\frac{\pi }{2}\) and \(\frac{3\pi }{2}\) angle rotated and mirror-symmetrical images.

The places of data points are also in accordance to that expected, since the randomized group (green dots) attaches \(\mathcal{H}\!\sim \!1\), \(\mathcal{C}_{JS}\!\sim \!0\) and \(\mathcal{F}\!\sim \!0\), the ordered (blue dots) has the minimum \(\mathcal H\) value, a low \(\mathcal{C}_{JS}\) value and a higher \(\mathcal{F}_{Lehmer}\) value. The multifractal group (red dots) attaches intermediate values of \(\mathcal{H}\) and \(\mathcal{C}_{JS}\) (being the highest of the three groups), and a high \(\mathcal{F}_{Lehmer}\) value similar to the ordered group, which shows the ability of the method to properly capture the rich structuring of these type of images.

4.2 Brownian surfaces

In contrast to the deterministic process generating an SMM, the Brownian surfaces are the product of a random process, so the resulting images are closer to real-world ones. They have been introduced by B. Mandelbrot [25] for simulating natural textures, and have since been successfully used to design computer-generated landscapes.

Here, a random variable, X(xy), is considered as the height of a surface at each point \((x,y)\!\in \!E\). For \(0\!<\!H\!<\!1\), an index-H fractional Brownian function \(X\!:\!{\mathbb {R}}^2\!\rightarrow \!{\mathbb {R}}\) is a Gaussian random function such that [11]: (i) with probability 1, \(X(0,0)\!=\!0\) and X(xy) is a continuous function of (xy); (ii) for (xy), \((h,k)\!\in \!{\mathbb {R}}^2\), the height increments \(X(x+h,y+k)\!-\!X(x,y)\) are normally distributed, with mean 0 and variance \((h^2+k^2)^H\!=\!||(h,k)||^{2H}\).

The set \(\big \{\big ((x,y),X(x,y)\big )\!:(x,y)\!\in \!{\mathbb {R}}^2\big \}\) is called an index-H fractional Brownian surface (H-fBs) and it is known that, with probability 1, it has both, Hausdorff and box-counting fractal dimension equal to \(3\!-\!H\) [11]. Thus, the closer the H value is to 1, the lower is its fractal dimension and so its level of roughness. On the contrary, the closer the H value is to 0, the higher the fractal dimension and the level of roughness. The index H is also known in the literature as the Hurst exponent and is indeed used as a roughness parameter. Three numerical simulated examples of H-fBs of size \(2^{10}\!\times \!2^{10}\) are illustrated in Figs. 8, 9 and 10, for \(H\!=0.9,0.5\) and 0.1, respectively.

Fig. 8
figure 8

a Numerical simulation of a H-fBs on the unit square for \(H\!=\!0.9\), b 3-dimensional rendering of (a)

Fig. 9
figure 9

a Numerical simulation of a H-fBs on the unit square for \(H\!=\!0.5\), b 3-dimensional rendering of (a)

Fig. 10
figure 10

a Numerical simulation of a H-fBs on the unit square for \(H\!=\!0.1\), b 3-dimensional rendering of (a)

Fig. 11
figure 11

CECP (a) and FECP (b) for 30 independent numerical simulations of H-fBs for each H, \(H\!\in \!\{0.1,0.2,\ldots ,0.9\}\), for \(D\!=\!5\) and \(\tau \!=\!1\). Mean and standard deviation for each set of realizations are also plotted as error bars

In this work, the information quantifiers involved were calculated for 30 independent simulations for each value of H, \(H\!\in \!\{0.1,0.2,\ldots ,0.9\}\). The results are placed in the CECP and FECP in Fig. 11(a) and (b) respectively, for the parameter values \(D\!=\!5\) and \(\tau \!=\!1\), together with the mean and standard deviations of each set of realizations. It can be appreciated that both planes discriminate well the groups taking into account the random component that all images have. The data points in the CECP fall in a curve according to the expected values of \(\mathcal{H}\!\sim \!1\) and \(\mathcal{C}\!\sim \!0\) for negatively correlated surfaces (\(H\!<\!0.5\)) whose images appear more ‘disordered’, while for positive correlated (\(H\!>\!0.5\)) the entropy \(\mathcal{H}\) decreases and the complexity \(\mathcal{C}_{JS}\) increases as H scales to 0.9. In the FECP instead, the data points locate in a straight line revealing what appears to be an inverse linear relationship between \(\mathcal H\) and \(\mathcal F\). Note that \(\mathcal F\) can reach values higher than \(\mathcal C\) ones for \(H\!=\!0.8,0.9\) and, on the contrary, lower values for \(H\!=\!0.1, 0.2\). This is because, as a local quantifier, \(\mathcal F\) is more sensitive to capture local correlations just characterized by the H exponent. This illustrates the relevance of adding the FECP to the set of tools for studying image complexity.

Fig. 12
figure 12

CECP (a) and FECP (b) for a single simulation for each H value in Fig. 11 along with the corresponding rotated (\(\frac{\pi }{2},\pi\) and \(\frac{3\pi }{2}\)) and mirror-symmetrical images, \(D\!=\!5\), \(\tau \!=\!1\)

Figures 12(a) and (b) depict the CECP and FECP respectively, for a single realization of each H value in Fig. 11 along with their rotations at \(\frac{\pi }{2}, \pi\) and \(\frac{3\pi }{2}\), and the mirror-symmetrical image. The invariance of the results against those transformations is here also outright.

Additionally, we register the execution times of the proposed method and the one by [48]. For both methods, the experiment reads 180 simulated images, with the Hurst exponent ranging from 0.1 to 0.9, and then computes the information theory quantifiers. The experiments were performed using MATLAB (R2023a) on a PC with Intel(R) Core(TM) i7-8700 CPU @3.2GHz, 6 cores, and 16 GB RAM. Table 1 shows the results. Our method appears less costly in terms of machine time for \(D>4\). The overall evaluation is that, for \(D=6\), reading the images and computing their information theory quantifiers takes approximately one minute and a half. Even for a large pattern length (\(D=8\)) it takes approximately 27 min.

Table 1 Execution time (in seconds) for reading 120 simulated images and computing the information-theory quantifiers for comparable pattern lengths, using the proposed method and a benchmark method described by [48]

4.3 Scale invariance

The delay parameter \(\tau\) can be regarded as a scale resolution factor when analyzing images. Figures 13, 14, and 15, depict the variation of \({\mathcal {H}}\), \({\mathcal {C}}_{JS}\), and \({\mathcal {F}}_{Lehmer}\), respectively, as a function of \(\tau\) for a fixed embedding dimension \(D\!=\!6\).

Initially, the range of quantifier values is larger for lower values of \(\tau\). This is expected, as high-resolution images are typically more easily distinguishable, and neighboring pixels exhibit minimal differences.

However, not all Brownian surfaces exhibit identical behavior. The quantifier values for standard Brownian motion (\(H\!=\!0.5\)) and anti-persistent fractional Brownian motion (\(H\!<\!0.5\)) are only marginally affected by the time delay for \(\tau \!>\! 2\). Conversely, the impact of the time delay is more pronounced for persistent processes.

Furthermore, each Hurst exponent is characterized by a specific value of \({\mathcal {H}}\), \({\mathcal {C}}_{JS}\), and \({\mathcal {F}}_{Lehmer}\) for each \(\tau\). In other words, the evolution of the quantifiers does not overlap. Consequently, despite the narrowing range of quantifier values as \(\tau\) increases, the characterization of Brownian surfaces remains valid for different \(\tau\) values. This renders the proposed method invariant to image resolution.

A closer examination of the figures reveals that all curves exhibit a more or less pronounced peak at \(\tau\) values that are multiples of 4. This peculiar behavior is intrinsically linked to the inner structure of the Hilbert curve, as every 4 steps in its path, it turns towards a previous neighbor (see Fig. 1). As \(\tau\) is a parameter related to the autocorrelation of the time series, these peaks are reasonably more pronounced for curves associated with higher values of H, which correspond to positively autocorrelated data.

Fig. 13
figure 13

Shannon entropy \(\mathcal H\) as a function of \(\tau\), for Hurst exponents \(H\!\in \!\{0.1,0.2,\ldots ,0.9\}\) and fixed embedding dimension \(D\!=\!6\)

Fig. 14
figure 14

Complexity \(\mathcal{C}_{JS}\) as a function of \(\tau\), for Hurst exponents \(H\!\in \!\{0.1,0.2,\ldots ,0.9\}\) and fixed embedding dimension \(D\!=\!6\)

Fig. 15
figure 15

Fisher Information Measure \(\mathcal{F}_{Lehmer}\) as a function of \(\tau\), for Hurst exponents \(H\!\in \!\{0.1,0.2,\ldots ,0.9\}\) and fixed embedding dimension \(D\!=\!6\)

4.4 Normalized Brodatz

We apply the method to a real-world surface dataset by taking images from the Normalized Brodatz Texture database (NBT) [1].Footnote 1 This database comprises 112 photos from the album ‘Textures: A Photographic Album for Artists and Designers’ [7]. The original Brodatz images have different background intensities and some of these texture images could be discriminated using only this feature (i.e., first-order statistics) without using texture information. Thus, they were enhanced by a normalization process which removed the greyscale background effect. Figure 16 shows 8 different images from the NBT database, selected as follow-up examples of the various techniques we implemented throughout the work. In addition, some of them coincide in having been analyzed in other works, such as [49], for a possible confrontation. The images of the NBT database are \(640\!\times \!640\) pixels in size and use 8 bits per pixel, allowing for 256 different shades of gray. They were cropped to a size of \(512\!\times \!512\) to be scanned by a Hilbert curve of level 9.

Fig. 16
figure 16

Eight Normalized Brodatz Textures selected from NBT database with their corresponding labels

Fig. 17
figure 17

CECP (a) and FECP (b) for the Normalized Brodatz set with \(D\!=\!8\) and \(\tau \! =\! 1\). The data points and labels of the eight selected textures (Fig. 16) are highlighted in red

Both planar representations, CECP and FECP, were calculated for the 112 images of the database. Figure 17 illustrates the locations in both planes for the case \(D\!=\!8\) and \(\tau \!=\!1\), where the data points corresponding to the 8 selected images are highlighted in red. It can be seen that they spread over the suggested representation planes, but the FECP discriminates the images better. Thus, it’s crucial to calculate both quantifiers. Regarding the selected textures in the FECP, for example, one can distinguish in Fig. 17, three groups according to their FIM values: D16, D71 and D93 (the lowest), D15 and D49 (the medium ones) and D44, D101 and D102 (the highest). Looking at these groups in Fig. 16, it is possible to relate them with different types of patterns: D44, D101 and D102 have a regular design with certain symmetry around a central point. Besides D101 and D102 display a periodic motif and so they have a higher value of entropy \(\mathcal H\) than D44. D15 and D49 seems to have a preferential direction of plot drawing, vertical and horizontal, respectively, being sharper in D49, and that is why its entropy is greater than that of D15. Finally, D71 and D93 are more ‘disordered’ designs and D16 has almost no contrasts being therefore close to the uniform distribution.

Fig. 18
figure 18

Entropy-Complexity-Fisher space of Normalized Brodatz Texture database. Embedding dimensions \(D_x\!\times \!D_y\!=\!2\!\times \!4\) and delays \(\tau _x\!=\!1\), \(\tau _y\!=\!1\) for Zunino & Ribeiro method ([48]) (red points); \(D\!=\!8\) and \(\tau \!=\!1\) for the proposed method (blue points)

In order to compare our results with those by Zunino & Rivero [48], we propose a three-dimensional representation of the information quantifiers. For each P the triple \(\big (\mathcal{H}[P],\mathcal{C}_{JS}[P],\mathcal{F}[P]\big )\) is located in the \(\mathcal{H}\!\times \!\mathcal{C}\!\times \!\mathcal{F}\) space. We also compute the FIM for the 112 pictures of the NBD, following the row-wise order used in [48] for scanning the images when computing the statistical complexity, by using here submatrices of embedding dimensions \(D_x\!\times \!D_y\!=\!2\!\times \!4\) and delays \(\tau _x\!=\!1\), \(\tau _y\!=\!1\). For that aim we used an embedding dimension \(D\!=\!8\) and \(\tau \!=\!1\) when applying the Hilbert curve method (our proposal). The results are depicted in Fig. 18. We observe that our approach reaches higher values for the permutation entropy and begins with lower values of statistical complexity and FIM. This figure shows that adding the FIM is useful to enhance image discriminationFootnote 2

4.5 Colored Brodatz

Fig. 19
figure 19

Colored Brodatz Textures selected from CBT database with their corresponding labels. They are colored versions of the textures in Fig. 16

In [1], it is also offered a Colored Brodatz Texture database (CBT), which consists of a colored version of the 112 Brodatz images, and has the advantage of retaining the rich textural content of the originals while having a wide variety of colors. Figure 19 shows the colored versions of the 8 selected pictures in Fig. 16. We perform the proposed method to these new images and the results are depicted in Fig. 20. Comparing the latter with Fig. 16 it can be seen that the location of the normalized and colored textures in the CECP and FECP representations do not differ significantly, as highlighted with the 8 pictures selected for this comparison. The only exception in the color version is one of the images (labeled D111 in both databases), which behaves as an outlier, probably due to the over-coloring effect. Therefore, the proposed method could be applied either to black and white or to color images.

Fig. 20
figure 20

Colored Brodatz set with with \(D\!=\!8\), \(\tau \!=\!1\) (a) CECP (b) EFCP

4.6 Rotated Brodatz

It is worthwhile to test the performance of the proposed method on textures rotated by angles other than just multiples of \(\frac{\pi }{2}\). The Signal and Image Processing Institute (SIPI) of the University of Southern CaliforniaFootnote 3 provides a set of images containing Brodatz textures rotated by seven different angles: 0, 30, 60, 90, 120, 150 and 200 degrees. For this evaluation, we selected eight different textures and their corresponding rotations: Raffia, Straw, Weave, Wool, Leather, Water, Pigskin, and Brick. This selection encompasses a variety of textures, making it representative of the entire Brodatz database. Figure 21 shows an example of them (Brick).

Fig. 21
figure 21

A Normalized Brodatz texture and its rotated versions. Each panel is labeled with the degree of rotation

Figure 22 illustrates the CECP and the FECP representations for \(D\!=\!8\) and \(\tau \!=\!1\) of the 8 selected textures and their corresponding rotations. It can be observed that the locations of data points in both planes are fairly invariant to rotations, considering the loss of information. That is to say, the discrepancy between some points within the same group is greater than in previous examples. This is because rotations by angles that are multiples of \(90^{\circ }\) map squares onto squares, and thus ‘all’ the information of the original image is necessarily preserved. In this case, an image is rotated and then both the original and rotated versions are cropped to the same square size. Therefore, while they do not share exactly the same information, they retain most of it.

Fig. 22
figure 22

CECP (a) and FECP (b) for Rotated Brodatz set with \(D\!=\!8\), \(\tau \!=\!1\)

5 Conclusions

We propose a novel method, based on information quantifiers, to distinguish between textures in images, leveraging the locality-preserving properties of the Hilbert curve as a means of reading them. One of the advantages of reading images using the Hilbert curve is that it avoids the directional biases found in conventional methods. Additionally, we introduce the calculation of the FIM (as a local quantifier) and the FECP representation to the existing literature, as a complementary indicator for analyzing complexity in two-dimensional patterns.

Briefly, the method is composed by three steps:

  • Read the image using the Hilbert curve.

  • Compute the probability density function according to the Bandt and Pompe symbolization method for embedding dimensions \(5\!\le \!D\!\le \!8\) and delay \(\tau \!=\!1\).

  • Compute the information theory quantifiers.

This method has been applied to simulated and real images with satisfactory results. Regarding the image of the SMM and its both ordered and randomized versions, the proposed method result efficient to distinguish between regular, random, and more complex designs. It results also a robust method against rotations and mirror-symmetry transformations, as shown in Fig. 7. Regarding the simulated H-fBs, we generated surfaces with different values of Hurst exponents, which occupied different locations in the CECP and FECP, and that also resulted invariant to rotation, mirror symmetry, and image scaling.

Regarding the real-world images (Normalized and colored Brodatz data-bases), we observe a similar location in the planes, indicating that the proposed method is robust to the different color settings. Finally, we provide evidence that our discrimination method is also invariant to rotations.

The invariance to rotation and color settings are important properties of an image discrimination method. Another advantage of our method is its natural extension to higher dimensions employing the n-dimensional Hilbert curve, as well as its possible adaptation to non-square images. Thus, the results presented in this paper will be useful for scientists working with diverse kinds of images and textures.

Future work includes the application of this method to the analysis of medical images, aimed at discriminating different pathologies. Also, the consideration of alternative definitions of entropy, like logical entropy, defined by Ellerman [9, 26], and its extensions to fuzzy logic [6], will be a promising line of research. Another extension of this project could include the use of alternative information theory quantifiers such as the Renyi or the Tsallis entropies.