Texture Analysis

Texture analysis describes an image region by the statistics of its local structure rather than by the layout of its contents. A patch of sand, woven cloth, bark or skin is recognized by what kinds of small structures it contains, how often they occur, and how they relate to one another at short distances, not by where any single grain or thread lies. The field’s methods compute such statistics, compare them, and use them to classify materials, to segment an image into regions of uniform texture, and to synthesize new samples of a texture. The main families are co-occurrence statistics of grey levels, energies of a bank of oriented band-pass filters, histograms of local binary codes, histograms of learned prototypes of filter responses (textons), and the responses of convolutional networks pooled without regard to position (Cimpoi, Maji and Vedaldi, 2015). This article explains what a texture statistic is, what the families measure, and, in a worked example on textures with known generating parameters, what each can and cannot tell apart.

Definition

There is no single formal definition of texture. Portilla and Simoncelli call the definition “rather imprecise” on their laboratory’s texture page and give the one most researchers agree on: an image of visual texture is spatially homogeneous and usually contains repeated structures, often with random variation in position, orientation or colour. Their 2000 paper makes this precise in the framework that, they write, Julesz’s work established: a texture is modelled as a homogeneous (stationary) random field on the pixel lattice, and a texture representation is a set of statistics, functions whose spatial averages over an image estimate expectations under that field (Portilla and Simoncelli, 2000, §1).

Texture analysis is then the computation and use of such statistics. A texture descriptor is the vector of statistics computed over a region; the region must be large enough for the averages to be stable but small enough to contain one texture.

Intuition

Two pictures of the same gravel never agree pixel by pixel, yet they look alike. What they share is a distribution: the grey levels that occur, how often a bright pixel has a bright neighbour two pixels to the right, how much of the image’s variation lies at fine scales and how much in vertical rather than horizontal structure. A texture statistic averages some local measurement over the whole region, so it is insensitive to where things are and sensitive to what kinds of things are there. This is the same pooling that turns local descriptors into a bag of visual words, and Cimpoi, Maji and Vedaldi (2015) describe it as “orderless pooling”: texture is traditionally described by pooling filter-bank responses without regard to their arrangement, since overall shape is usually unimportant for a texture.

The families differ in which local measurement they average:

• First-order statistics: the image histogram of grey levels, or its moments. They ignore spatial structure entirely.
• Second-order statistics of pixel pairs: the co-occurrence matrix, which counts how often grey level $i$ appears at a fixed offset from grey level $j$, and the autocorrelation or, equivalently, the power spectrum.
• Filter-bank energies: the image is filtered with band-pass filters at several scales and orientations, such as Gabor filters or derivatives of Gaussians, and the squared responses are averaged. Each filter measures how much of the image’s variation lies in its frequency band.
• Codes and prototypes: each pixel’s neighbourhood is mapped to one of a finite set of labels, by comparisons with the centre pixel (local binary patterns) or by assigning the vector of filter responses to its nearest cluster centre (textons), and the region is described by the histogram of labels.

The central question is which statistics are enough. Julesz posed it for human perception: are two textures that agree in all statistics up to some order indistinguishable at a glance? The question frames the worked example below, which builds pairs of textures that share their autocorrelation and power spectrum but differ in their grey-level distribution and in the joint distribution of pixel pairs.

Formal Definition

Texture as a stationary random field

Let $X$ be a homogeneous random field on the lattice. Portilla and Simoncelli (2000, eq. 1) state the hypothesis they attribute to Julesz (1962), in the form later authors reformulated it: there is a set of constraint functions $\phi_k$ such that two fields $X$ and $Y$ with $E[\phi_k(X)] = E[\phi_k(Y)]$ for every $k$ produce samples that look the same. A texture representation is a choice of $\phi_k$, and an image’s descriptor is the vector of spatial averages that estimate them, which requires that the averages over one finite image be close to the expectations; Portilla and Simoncelli call this property practical ergodicity (§1.1).

Grey-level co-occurrence matrices

Quantize the image $I$ to $L$ grey levels, $q(\mathbf{x}) \in \{0, \ldots, L-1\}$, and fix an offset $\mathbf{d} = (\Delta r, \Delta c)$. The grey-level co-occurrence matrix (GLCM) counts the pixel pairs at that offset,

$$P_{\mathbf{d}}(i, j) = \big|\{\mathbf{x} : q(\mathbf{x}) = i,\; q(\mathbf{x} + \mathbf{d}) = j\}\big| ,$$

over all pairs whose both pixels lie in the region, and is normalized to sum to 1. It is the empirical joint distribution of a pixel and its neighbour at offset $\mathbf{d}$. Summary features, which scikit-image’s graycoprops calls properties, reduce each matrix to a number:

$$\text{contrast} = \sum_{i,j} P(i,j)\,(i-j)^2, \qquad \text{homogeneity} = \sum_{i,j} \frac{P(i,j)}{1+(i-j)^2}, \qquad \text{energy} = \Big(\sum_{i,j} P(i,j)^2\Big)^{1/2},$$ $$\text{correlation} = \sum_{i,j} P(i,j)\,\frac{(i-\mu_i)(j-\mu_j)}{\sigma_i \sigma_j},$$

with $\mu_i, \sigma_i$ and $\mu_j, \sigma_j$ the means and standard deviations of the two marginals of $P$. scikit-image’s documentation of graycomatrix cites Haralick, Shanmugam and Dinstein (1973) and a tutorial by Hall-Beyer; the formulas above are the ones in scikit-image’s source, which also offers dissimilarity, the angular second moment (energy squared), entropy, and the mean, variance and standard deviation of the matrix. Contrast and homogeneity measure how far apart the paired grey levels tend to be, energy how concentrated the matrix is, and correlation the linear dependence between a pixel and its neighbour. A descriptor can combine several offsets, for example a few distances in four directions, as the example below does.

The autocorrelation and the power spectrum

For a stationary field with zero mean and unit variance, the autocorrelation is $\rho(\mathbf{d}) = E[X(\mathbf{x})X(\mathbf{x}+\mathbf{d})]$, and its Fourier transform is the power spectrum $S(\mathbf{f})$, the distribution of the field’s variance over spatial frequencies (see Fourier transform, whose list of uses includes summarizing a texture by the energy in rings and wedges of the power spectrum). The GLCM contains more than $\rho$: it is the whole joint distribution of the pair, not only its correlation. For a Gaussian field, however, the joint distribution of any pair is bivariate normal with correlation $\rho(\mathbf{d})$, so the GLCM of a quantized Gaussian field is a function of $\rho(\mathbf{d})$ alone, and the example computes it exactly from $\rho$.

Filter-bank energies and the power spectrum

The Gabor filter article covers the filter itself: its frequency and orientation tuning, quadrature pairs, and how a bank spaces its scales and orientations. Let $h_k$ be the filters of a bank, with transfer functions $H_k(\mathbf{f})$, and let $r_k = h_k \ast I$. The energy feature of filter $k$ is the mean of $|r_k|^2$ over the region. The following is our derivation, for an $N \times N$ periodic image and the discrete Fourier transform without normalization (NumPy’s convention). By Parseval’s theorem and the convolution theorem,

$$\frac{1}{N^2}\sum_{\mathbf{x}} |r_k(\mathbf{x})|^2 = \frac{1}{N^4} \sum_{\mathbf{f}} |H_k(\mathbf{f})|^2 \, |\hat I(\mathbf{f})|^2 ,$$

a linear function of the image’s periodogram $|\hat I|^2$. Taking expectations over a stationary field whose expected periodogram is $N^2 S(\mathbf{f})$ (so that its variance is $N^{-2}\sum_{\mathbf{f}} S(\mathbf{f})$),

$$E\Big[\frac{1}{N^2}\sum_{\mathbf{x}} |r_k(\mathbf{x})|^2\Big] = \frac{1}{N^2}\sum_{\mathbf{f}} |H_k(\mathbf{f})|^2\, S(\mathbf{f}) .$$

The expected mean energy of every filter is therefore a weighted integral of the power spectrum, whatever else the field’s statistics are. Two textures with the same power spectrum have the same expected mean energies in every bank, and no classifier of mean energies alone can tell them apart better than the sampling noise allows. Statistics of the energy beyond its mean, such as its variance over the region, involve fourth-order moments of the image and can. A further consequence, also our derivation: if the field is Gaussian and the filter is complex with negligible response at negative frequencies, as a complex Gabor filter tuned away from zero frequency is, $r_k$ is a circular complex Gaussian and $|r_k|^2$ is exponentially distributed, so its standard deviation equals its mean and its coefficient of variation is 1.

Textons

Leung and Malik (2001, §2) filter each image with a bank of 48 filters, 36 elongated filters at 6 orientations, 3 scales and 2 phases, 8 centre-surround filters (difference-of-Gaussian filters in the text, derivative filters in the caption of Fig. 4) and 4 low-pass Gaussians, so that each pixel becomes a 48-dimensional vector, and cluster these vectors with k-means. The cluster centres are the textons; each pixel is labelled with its nearest texton, and, in their material recognition with 3D textons (§5), a material is described by the histogram of labels, compared with a $\chi^2$ statistic. Malik, Belongie, Leung and Shi (2001, §1.2) present this as an operational definition of the textons that Julesz (1981) had proposed as the units of preattentive texture perception and had described only qualitatively, as oriented segments, crossings and terminators in binary line patterns. The texton histogram is a bag of visual words whose local descriptor is a filter-response vector.

Local binary patterns

A local binary pattern codes each pixel by which of its neighbours on a circle are at least as bright as the centre, and a region by the histogram of codes, optionally after mapping codes to rotation-invariant “uniform” labels. The local binary patterns article covers the operator, its mappings and its invariances: it is unchanged by any strictly increasing grey-level map when its neighbours lie on pixel centres, and it discards contrast.

Properties

What a family can distinguish follows from what it measures. First-order statistics cannot separate textures with the same grey-level distribution, however different their spatial structure. Mean filter-bank energies cannot separate textures with the same power spectrum (derived above). The GLCM at a set of offsets cannot separate textures whose pixel pairs at those offsets have the same joint distributions. Portilla and Simoncelli (2000, §2.2) note that correlation measurements alone had long been known to miss the structure of many natural textures; their model adds marginal statistics of the pixels, correlations of the magnitudes of complex wavelet coefficients across neighbouring positions, orientations and scales, and relative phases across scales, and they show for each group of statistics a texture whose synthesis fails without it.

The Julesz conjecture. Two sources date its history differently. Portilla and Simoncelli (2000) write that Julesz (1962) hypothesized that $N$th-order joint statistics of the pixels, for some unspecified $N$, partition textures into classes that cannot be told apart preattentively; that Julesz et al. (1973) later proposed that second-order ($N = 2$) statistics suffice; and that this was disproved by pairs of textures with identical second- and even third-order statistics that look different (Caelli and Julesz, 1978; Julesz et al., 1978). The laboratory’s texture page instead names Julesz (1962) as the origin of the conjecture that humans cannot tell apart textures with identical second-order statistics, and says that he disproved it in a 1973 paper with counterexamples. Both agree that the second-order form of the conjecture fails for human vision.

Invariance. A descriptor for materials may need to ignore illumination, contrast, rotation and scale. Each family achieves this differently: GLCM features can be averaged or pooled over directions, filter-bank features can be pooled over orientations or normalized by the total energy, LBP codes are invariant to monotonic grey-level changes and have rotation-invariant mappings. Leung and Malik (2001, §1) point out a limit that brightness normalization and intensity transforms such as histogram matching do not remove: for rough surfaces, a ridge that appears as a light–dark transition under one illumination appears as a dark–light transition under another, so changes of lighting and viewpoint alter the texture itself; their 3D textons model a material’s appearance across several lighting and viewing directions at once.

Scale and window size. A descriptor estimates expectations from one finite region, so it fluctuates from region to region, and more for textures whose correlation length is large compared with the region (fewer independent samples). For segmentation the window has a further cost: Malik et al. (2001, §1) show that a window straddling the boundary between two textures contains large responses of filters aligned with the edge, so a segmentation by histogram distances tends to carve the boundary out as a thin region of its own.

Quantization and offsets in the GLCM. The GLCM has $L^2$ entries per offset, so images are quantized to a small number of levels; scikit-image’s docstring suggests binning the image rather than using a large number of levels. In scikit-image 0.25.2, the offset for distance $d$ and angle $\theta$ is computed in _glcm_loop as $(\Delta r, \Delta c) = (\operatorname{round}(d\sin\theta), \operatorname{round}(d\cos\theta))$, so $\theta = \pi/2$ pairs each pixel with the pixel $d$ rows below it. The docstring describes $\pi/2$ as the pixel above; its own example output for $\pi/2$ matches the code, not the description. The rounding also means that a diagonal offset of distance 2 is $(1, 1)$, the same pair as distance 1, so the GLCM features for distances 1 and 2 at $45°$ and $135°$ are identical, and distance 4 on a diagonal is $(3, 3)$, a separation of $4.24$ pixels.

Examples

The example builds textures whose statistics are known exactly and asks how each family’s features track the generating parameters and what each can separate. Each texture is a periodic $128 \times 128$ field made by filtering white noise with an anisotropic Gaussian kernel: standard deviation $\ell\sqrt{a}$ along a direction $\varphi$ and $\ell/\sqrt{a}$ across it, where $\ell$ is the correlation length, $a$ the anisotropy and $\varphi$ the orientation, measured counter-clockwise from the column axis with rows pointing down. The noise is of two kinds. Gaussian white noise gives a Gaussian random field. Shot noise, independent Poisson counts of 0.02 spots per pixel minus their mean and scaled to unit variance, gives a field of randomly placed blurred spots with exactly the same autocorrelation and power spectrum, since both noises are white with unit variance, but with a skewed distribution and non-Gaussian pair statistics. Both kinds are scaled to unit variance. The exact model gives the autocorrelation from the kernel’s transfer function $K$ (as $|K|^2$), the shot noise’s skewness from the discrete kernel and the Poisson third cumulant, the GLCM of a quantized Gaussian pair by numerical integration of the bivariate normal, and each filter’s expected mean energy from the formula derived above.

Four families describe each image: a 16-bin grey-level histogram; GLCM properties (contrast, homogeneity, energy, correlation) from scikit-image’s graycomatrix and graycoprops on 16 levels, equal-probability for a standard normal distribution, at distances 1, 2 and 4 in four directions; the log mean energy of 18 complex Gabor filters from scikit-image’s gabor_kernel (frequencies 0.25, 0.125 and 0.0625 cycles per pixel, six orientations, one-octave bandwidth), applied by circular convolution since the fields are periodic; and LBP histograms with the riu2 mapping at $(P, R) = (8, 1)$ and $(16, 2)$ on the field mapped to 8 bits. The GLCM and Gabor features appear in two forms: directional, as computed, and sorted, where the values for the different directions of each property and scale are sorted, which makes them unchanged by any permutation of the directions (our choice of pooling). A third Gabor form adds each filter’s coefficient of variation (CV) of energy over the image. Part A checks the hand-written GLCM against scikit-image, the FFT filtering against scipy.ndimage.convolve with wrap-around, and the generated data against the exact model. The script ran in about 14 seconds on one CPU thread with NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2 and scikit-learn 1.7.2; it takes the seed as an argument and the output below is seed 0, run twice with identical output.

import os, sys
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[v] = "1"                                   # one BLAS thread, before NumPy loads
import numpy as np
from scipy import ndimage, stats
from skimage.feature import graycomatrix, graycoprops, local_binary_pattern
from skimage.filters import gabor_kernel
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

SEED = int(sys.argv[1]) if len(sys.argv) > 1 else 0      # base seed of every random field
N, LAM, L = 128, 0.02, 16                                 # image side, shot-noise density, GLCM levels
FY, FX = np.meshgrid(np.fft.fftfreq(N), np.fft.fftfreq(N), indexing="ij")  # cycles/pixel (rows, cols)

# ---------- textures with known generating parameters ----------
def transfer(ell, a, phi):
    """Transfer function of an anisotropic Gaussian kernel: standard deviation ell*sqrt(a) along
    the direction phi (counter-clockwise from the column axis, rows pointing down) and ell/sqrt(a)
    across it; scaled so that filtered unit-variance white noise has unit variance."""
    s1, s2 = ell * np.sqrt(a), ell / np.sqrt(a)
    u = FX * np.cos(phi) - FY * np.sin(phi)               # frequency component along phi
    v = FX * np.sin(phi) + FY * np.cos(phi)               # ... and across it
    K = np.exp(-2 * np.pi**2 * (s1**2 * u**2 + s2**2 * v**2))
    return K * np.sqrt(N * N / np.sum(K**2))

def field(kind, ell, a, phi, rng):
    """Periodic N x N texture: white noise (Gaussian, or centred Poisson 'shot' noise of density LAM
    spots per pixel, scaled to unit variance) filtered with the kernel. Both kinds have the same
    autocorrelation and power spectrum; only the Gaussian one is a Gaussian random field."""
    if kind == "gaussian":
        w = rng.standard_normal((N, N))
    else:
        w = (rng.poisson(LAM, (N, N)) - LAM) / np.sqrt(LAM)
    return np.real(np.fft.ifft2(np.fft.fft2(w) * transfer(ell, a, phi)))

def rho_exact(ell, a, phi):
    """Exact autocorrelation of either kind at every (row, col) offset (periodic)."""
    c = np.real(np.fft.ifft2(transfer(ell, a, phi) ** 2))
    return c / c[0, 0]

def skew_exact(ell, a, phi):
    """Exact skewness of the shot-noise field (third cumulant of Poisson counts = LAM)."""
    k = np.real(np.fft.ifft2(transfer(ell, a, phi)))
    return np.sum(k**3) / np.sum(k**2) ** 1.5 / np.sqrt(LAM)

# ---------- descriptor families ----------
EDGES = stats.norm.ppf(np.arange(1, L) / L)               # 16 equal-probability levels of N(0, 1)
quantize = lambda z: np.searchsorted(EDGES, z).astype(np.uint8)
DIST, ANG = [1, 2, 4], [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4]
OFFS = [[(int(round(np.sin(t) * d)), int(round(np.cos(t) * d))) for t in ANG] for d in DIST]
PROPS = ["contrast", "homogeneity", "energy", "correlation"]

def glcm_hand(q, dr, dc):
    """Co-occurrence counts P[i, j]: level i at (r, c), level j at (r + dr, c + dc), pairs inside the image."""
    H, W = q.shape
    r0, r1, c0, c1 = max(0, -dr), min(H, H - dr), max(0, -dc), min(W, W - dc)
    a, b = q[r0:r1, c0:c1].astype(int), q[r0 + dr:r1 + dr, c0 + dc:c1 + dc].astype(int)
    return np.bincount((a * L + b).ravel(), minlength=L * L).reshape(L, L)

def props_hand(P):
    P = P / P.sum()
    i, j = np.ogrid[0:L, 0:L]
    mi, mj = np.sum(i * P), np.sum(j * P)
    si, sj = np.sqrt(np.sum(P * (i - mi) ** 2)), np.sqrt(np.sum(P * (j - mj) ** 2))
    return {"contrast": np.sum(P * (i - j) ** 2), "homogeneity": np.sum(P / (1 + (i - j) ** 2)),
            "energy": np.sqrt(np.sum(P**2)), "correlation": np.sum(P * (i - mi) * (j - mj)) / (si * sj)}

def glcm_features(q):
    """(4 properties, 3 distances, 4 angles) from scikit-image."""
    P = graycomatrix(q, DIST, ANG, levels=L, normed=True)
    return np.array([graycoprops(P, p) for p in PROPS])

FREQS, THETAS = [0.25, 0.125, 0.0625], np.arange(6) * np.pi / 6
BANK = [gabor_kernel(f, theta=t) for f in FREQS for t in THETAS]   # bandwidth 1 octave (default)

def tf(k):
    """Transfer function of circular convolution with kernel k (centre at the origin)."""
    h = np.zeros((N, N), complex)
    h[:k.shape[0], :k.shape[1]] = k
    return np.fft.fft2(np.roll(h, (-(k.shape[0] // 2), -(k.shape[1] // 2)), axis=(0, 1)))
HB = [tf(k) for k in BANK]

def gabor_energy(z):
    """|response|^2 of the 18 complex Gabor filters (circular convolution; the fields are periodic)."""
    Z = np.fft.fft2(z)
    return np.array([np.abs(np.fft.ifft2(Z * H)) ** 2 for H in HB])

def energy_pred(ell, a, phi):
    """Expected mean energy of each filter: (1/N^2) sum_f |K(f)|^2 |H(f)|^2 (Parseval)."""
    K2 = transfer(ell, a, phi) ** 2
    return np.array([np.sum(K2 * np.abs(H) ** 2) / N**2 for H in HB])

def lbp_hist(z):
    """riu2 histograms for (P, R) = (8, 1) and (16, 2) on the field mapped to 8 bits, border dropped."""
    u8 = np.clip(np.round(128 + 40 * z), 0, 255).astype(np.uint8)
    h = []
    for P, R in ((8, 1), (16, 2)):
        c = local_binary_pattern(u8, P, R, "uniform")[3:-3, 3:-3].astype(int)
        h.append(np.bincount(c.ravel(), minlength=P + 2) / c.size)
    return np.concatenate(h)

def describe(z):
    q = quantize(z)
    g = glcm_features(q)                                  # (4, 3, 4)
    E = gabor_energy(z)
    m, cv = E.mean(axis=(1, 2)), E.std(axis=(1, 2)) / E.mean(axis=(1, 2))
    lm = np.log(m).reshape(3, 6)                          # (frequency, orientation)
    srt = lambda x: -np.sort(-x, axis=-1)                 # sort over directions: rotation-pooled
    return {"grey histogram (16)": np.bincount(q.ravel(), minlength=L) / q.size,
            "GLCM directional (48)": g.ravel(),
            "GLCM sorted (48)": srt(g).ravel(),
            "Gabor directional (18)": lm.ravel(),
            "Gabor sorted (18)": srt(lm).ravel(),
            "Gabor sorted + CV (36)": np.r_[srt(lm).ravel(), srt(cv.reshape(3, 6)).ravel()],
            "LBP riu2 8,1+16,2 (28)": lbp_hist(z)}

# ---------- exact model of the quantized Gaussian pair ----------
XG = np.linspace(-8, 8, 32001)
def pair_probs(rho):
    """P[i, j] of (quantize(X), quantize(Y)) for a standard bivariate normal with correlation rho."""
    e = np.r_[-np.inf, EDGES, np.inf]
    s = np.sqrt(1 - rho**2)
    cond = stats.norm.cdf((e[None, :] - rho * XG[:, None]) / s)          # P(Y < e_j | X = x)
    colp = np.diff(cond, axis=1) * stats.norm.pdf(XG)[:, None]          # density in x of Y-bin j
    idx = np.searchsorted(EDGES, XG)
    return np.array([np.trapezoid(colp * (idx == i)[:, None], XG, axis=0) for i in range(L)])

# ---------- A. checks: hand code against scikit-image and SciPy; data against the model ----------
print("A. checks")
rng = np.random.default_rng([SEED, 0])
dmax, pmax = 0, 0.0
for kind in ("gaussian", "shot"):
    q = quantize(field(kind, 3, 3, np.pi / 6, rng))
    Psk = graycomatrix(q, DIST, ANG, levels=L)
    for di in range(3):
        for ai in range(4):
            Ph = glcm_hand(q, *OFFS[di][ai])
            dmax = max(dmax, int(np.abs(Ph.astype(int) - Psk[:, :, di, ai].astype(int)).max()))
            ph = props_hand(Ph)
            pmax = max(pmax, max(abs(ph[p] - graycoprops(Psk[:, :, di:di + 1, ai:ai + 1], p)[0, 0]) for p in PROPS))
print(f"  hand GLCM vs graycomatrix: largest count difference {dmax}; hand properties vs graycoprops: {pmax:.1e}")
z = field("gaussian", 3, 3, np.pi / 6, rng)
E = gabor_energy(z)
gd = 0.0
for k, (kern, e) in enumerate(zip(BANK, E)):
    r = ndimage.convolve(z, kern.real, mode="wrap") + 1j * ndimage.convolve(z, kern.imag, mode="wrap")
    gd = max(gd, np.abs(np.abs(r) ** 2 - e).max())
print(f"  FFT Gabor energy vs scipy.ndimage.convolve (wrap): largest difference {gd:.1e}; "
      f"largest kernel {max(k.shape[0] for k in BANK)} px")
P0 = pair_probs(0.5)
print(f"  quantized-pair model: total probability {P0.sum():.6f}, largest marginal error {np.abs(P0.sum(1) - 1 / L).max():.1e}")
ell, a, phi = 3, 3, np.pi / 6
rho = rho_exact(ell, a, phi)
off = [(dr, dc) for dr in range(-8, 9) for dc in range(-8, 9)]
for kind in ("gaussian", "shot"):
    acs, mom = [], []
    for s in range(20):
        z = field(kind, ell, a, phi, np.random.default_rng([SEED, 1, s, kind == "shot"]))
        zc = z - z.mean()
        acs.append(np.real(np.fft.ifft2(np.abs(np.fft.fft2(zc)) ** 2)) / np.sum(zc**2))
        mom.append((z.mean(), z.var(), stats.skew(z.ravel())))
    ac = np.mean(acs, 0)
    err = max(abs(ac[dr, dc] - rho[dr, dc]) for dr, dc in off)
    m = np.mean(mom, 0)
    print(f"  {kind:8s} l=3 a=3 phi=30: |autocorr - exact| over |offsets| <= 8 (20 fields) {err:.3f}; "
          f"mean {m[0]:+.3f} var {m[1]:.3f} skewness {m[2]:+.3f} (exact {0 if kind == 'gaussian' else skew_exact(ell, a, phi):+.3f})")
print(f"  exact rho along phi at 1, 2, 4 px: "
      + ", ".join(f"{rho[int(round(-np.sin(phi) * d)), int(round(np.cos(phi) * d))]:.3f}" for d in (1, 2, 4))
      + "; across phi: " + ", ".join(f"{rho[int(round(-np.cos(phi) * d)), int(round(-np.sin(phi) * d))]:.3f}" for d in (1, 2, 4)))

Output:

A. checks
  hand GLCM vs graycomatrix: largest count difference 0; hand properties vs graycoprops: 2.2e-16
  FFT Gabor energy vs scipy.ndimage.convolve (wrap): largest difference 2.2e-15; largest kernel 55 px
  quantized-pair model: total probability 1.000000, largest marginal error 1.2e-04
  gaussian l=3 a=3 phi=30: |autocorr - exact| over |offsets| <= 8 (20 fields) 0.033; mean +0.014 var 0.967 skewness -0.036 (exact +0.000)
  shot     l=3 a=3 phi=30: |autocorr - exact| over |offsets| <= 8 (20 fields) 0.013; mean -0.007 var 0.984 skewness +0.872 (exact +0.887)
  exact rho along phi at 1, 2, 4 px: 0.973, 0.953, 0.883; across phi: 0.937, 0.660, 0.340

The checks pass to rounding error: the hand-counted co-occurrence matrices equal scikit-image’s at every offset, the FFT filtering equals spatial convolution, and the quantized-pair model integrates to one with each marginal within $1.2 \times 10^{-4}$ of $1/16$. The generated fields have the planned properties: over 20 fields the sample autocorrelation is within 0.033 (Gaussian) and 0.013 (shot) of the exact one at every offset up to 8 pixels, both kinds have mean near 0 and variance near 1, and the Gaussian fields are symmetric (skewness $-0.036$) while the shot-noise fields are skewed ($+0.872$ against an exact $+0.887$). The last line shows the anisotropy for $\ell = 3$, $a = 3$: four pixels along $\varphi$ the correlation is still 0.883, four pixels across it 0.340.

Part B measures, for nine combinations of the parameters and 40 fields of each kind per combination, how far each family’s averaged features lie from the exact model. The GLCM columns give the largest error of the correlation feature over the 12 offsets and the largest error of contrast relative to its exact value; the Gabor columns give the largest error of mean energy over the 18 filters, relative to the largest filter energy; CV is the coefficient of variation of energy averaged over the filters; the LBP columns give the share of uniform patterns at $(8, 1)$.

# ---------- B. how each family tracks the generating parameters ----------
CONDS = [(1.5, 1, 0), (3, 1, 0), (6, 1, 0), (1.5, 3, 0), (3, 3, 0), (6, 3, 0), (1.5, 3, 30), (3, 3, 30), (6, 3, 30)]
print("\nB. 40 fields per row; GLCM: largest |measured - exact| correlation, and |measured - exact| / exact contrast, over 3 distances x 4 angles;"
      "\n   Gabor: largest |measured - exact| mean energy over 18 filters, relative to the largest; CV of energy; LBP: uniform share")
print(f"  {'l':>4s}{'a':>3s}{'phi':>4s} | {'GLCM corr':>9s}{'contr.':>9s}{'shot corr':>10s}{'contr.':>9s} | "
      f"{'Gabor':>6s}{'shot':>6s}{'CV':>6s}{'shot CV':>8s} | {'LBP u':>6s}{'shot':>6s}")
for ci, (ell, a, phd) in enumerate(CONDS):
    phi = np.deg2rad(phd)
    rho = rho_exact(ell, a, phi)
    pred_c, pred_k = np.zeros((3, 4)), np.zeros((3, 4))
    for di in range(3):
        for ai in range(4):
            pp = props_hand(pair_probs(rho[OFFS[di][ai]]))
            pred_c[di, ai], pred_k[di, ai] = pp["correlation"], pp["contrast"]
    ep = energy_pred(ell, a, phi)
    row = []
    for kind in ("gaussian", "shot"):
        g, e, cvs, u = [], [], [], []
        for s in range(40):
            z = field(kind, ell, a, phi, np.random.default_rng([SEED, 2, ci, s, kind == "shot"]))
            g.append(glcm_features(quantize(z)))
            E = gabor_energy(z)
            e.append(E.mean(axis=(1, 2)))
            cvs.append(np.mean(E.std(axis=(1, 2)) / E.mean(axis=(1, 2))))
            u.append(1 - lbp_hist(z)[9])
        g = np.mean(g, 0)
        row.append((np.abs(g[3] - pred_c).max(), (np.abs(g[0] - pred_k) / pred_k).max(),
                    np.abs(np.mean(e, 0) - ep).max() / ep.max(), np.mean(cvs), np.mean(u)))
    (gc, gk, ge, gcv, gu), (sc, sk, se, scv, su) = row
    print(f"  {ell:4.1f}{a:3d}{phd:4d} | {gc:9.3f}{gk:9.3f}{sc:10.3f}{sk:9.3f} | {ge:6.3f}{se:6.3f}{gcv:6.2f}{scv:8.2f} | {gu:6.3f}{su:6.3f}")

Output:

B. 40 fields per row; GLCM: largest |measured - exact| correlation, and |measured - exact| / exact contrast, over 3 distances x 4 angles;
   Gabor: largest |measured - exact| mean energy over 18 filters, relative to the largest; CV of energy; LBP: uniform share
     l  a phi | GLCM corr   contr. shot corr   contr. |  Gabor  shot    CV shot CV |  LBP u  shot
   1.5  1   0 |     0.012    0.008     0.043    0.237 |  0.042 0.032  0.98    1.19 |  0.957 0.983
   3.0  1   0 |     0.003    0.008     0.003    0.035 |  0.027 0.041  0.98    1.10 |  0.987 0.990
   6.0  1   0 |     0.011    0.050     0.008    0.017 |  0.082 0.047  0.98    1.03 |  0.996 0.996
   1.5  3   0 |     0.006    0.003     0.055    0.238 |  0.011 0.015  0.98    1.18 |  0.958 0.977
   3.0  3   0 |     0.002    0.021     0.007    0.051 |  0.023 0.045  0.99    1.11 |  0.985 0.988
   6.0  3   0 |     0.007    0.013     0.012    0.021 |  0.004 0.016  0.97    1.02 |  0.993 0.993
   1.5  3  30 |     0.007    0.009     0.054    0.239 |  0.014 0.013  0.99    1.17 |  0.918 0.953
   3.0  3  30 |     0.020    0.022     0.007    0.026 |  0.029 0.070  0.98    1.09 |  0.973 0.976
   6.0  3  30 |     0.012    0.017     0.013    0.033 |  0.011 0.013  0.97    1.02 |  0.989 0.990

Co-occurrence statistics track the autocorrelation for Gaussian fields. For Gaussian fields the measured GLCM correlation is within 0.012 of the exact value in eight rows and 0.020 in the ninth ($\ell = 3$, $a = 3$, $\varphi = 30°$), and contrast within 5.0% of its exact value in every row, so, to within these errors, the GLCM features vary with $\ell$, $a$ and $\varphi$ as the autocorrelation at each offset does. The shot-noise fields have the same autocorrelation, yet at $\ell = 1.5$ their GLCM departs from the Gaussian model in all three groups: correlation by 0.043–0.055 and contrast by 23.7–23.9%. The pair distribution of sparse spots is not bivariate normal, and the GLCM, which records the whole pair distribution, sees it. The planned hypothesis was that the departure in contrast shrinks as $\ell$ grows in each group of $(a, \varphi)$, because more spots overlap within a correlation length and the field becomes more nearly Gaussian. It held for $a = 1$, $\varphi = 0$ (0.237, 0.035, 0.017) and for $a = 3$, $\varphi = 0$ (0.238, 0.051, 0.021). It failed for $a = 3$, $\varphi = 30°$, where the error falls from 0.239 to 0.026 and then rises to 0.033 at $\ell = 6$. In our reading, at $\ell = 6$ the shot-noise errors (0.017, 0.021, 0.033) are of the same size as the Gaussian ones in the same rows (0.050, 0.013, 0.017), which reflect only sampling scatter, so the step from $\ell = 3$ to $\ell = 6$ is too small for 40 fields to resolve; the large drop from $\ell = 1.5$ to $\ell = 3$ is clear in every group.

Mean filter energies follow the spectrum, for both kinds. The measured mean energies match the spectrum prediction to within 8.2% of the largest filter energy for Gaussian fields and 7.0% for shot noise, in every row: the two kinds, with the same spectrum, have the same expected energies, as derived. The coefficient of variation separates them: 0.97–0.99 for Gaussian fields, close to the value 1 derived above for an exponentially distributed energy, and 1.02–1.19 for shot noise, larger in every row and most at $\ell = 1.5$, where the spots are sparsest relative to their size.

LBP tracks smoothness. The share of uniform patterns rises with $\ell$ in every group for both kinds, as planned, from 0.918–0.958 at $\ell = 1.5$ to 0.989–0.996 at $\ell = 6$ for Gaussian fields: a smoother field has fewer sign alternations around a circle of radius 1. The LBP article measures the same share on independent noise and on images bundled with scikit-image. Of the planned checks in Part B, those on the Gaussian GLCM contrast (within 6% of the model), on the mean energies (within 10% for both kinds), on the coefficient of variation (within 0.05 of 1 for Gaussian fields, larger for shot noise) and on the uniform share held in every row. The check on the Gaussian GLCM correlation (within 0.02 of the model, limit included) failed in the row $\ell = 3$, $a = 3$, $\varphi = 30°$: its error prints as 0.020, and printing the same value with more digits gives 0.0201, just above the limit; it held in the other eight rows. The shot-noise departure at $\ell = 1.5$ held in all three groups, and its decrease with $\ell$ failed in one group, as described above.

Part C classifies. In C1, six classes of Gaussian fields differ in correlation length (2.5, 3, 3.6) and anisotropy (1, 1.5), and every image has a random orientation. In C2, the two classes are Gaussian and shot noise with the same spectrum, with $\ell$ and $a$ drawn at random for every image from $\{2, 3, 4.5\} \times \{1, 2.5\}$. Each descriptor is standardized with the training set’s statistics and classified by the nearest neighbour (1-NN) and by linear discriminant analysis (LDA), scikit-learn’s LinearDiscriminantAnalysis with its default deterministic solver. The C1 classes were moved closer once before the pilot, from $\ell \in \{2, 3, 4.5\}$ and $a \in \{1, 2.5\}$, because the first run of the planning script, on another seed, saturated; the shown classes were not changed further.

# ---------- C. classification by each family ----------
def classify(Xtr, ytr, Xte, yte):
    out = []
    for clf in (KNeighborsClassifier(1), LinearDiscriminantAnalysis()):
        sc = StandardScaler().fit(Xtr)
        out.append(np.mean(clf.fit(sc.transform(Xtr), ytr).predict(sc.transform(Xte)) == yte))
    return out

def experiment(name, classes, n):
    feats = {"train": [], "test": []}
    ys = {"train": [], "test": []}
    for split, sid in (("train", 0), ("test", 1)):
        for k, (kind, ells, aas) in enumerate(classes):
            for i in range(n):
                rng = np.random.default_rng([SEED, 3, name.startswith("C2"), sid, k, i])
                ell, a = rng.choice(ells), rng.choice(aas)
                phi = rng.uniform(0, np.pi)                   # random orientation for every image
                feats[split].append(describe(field(kind, ell, a, phi, rng)))
                ys[split].append(k)
    print(f"\n{name}: {len(classes)} classes, {n} training and {n} test images per class, random orientation")
    print(f"  {'descriptor':26s}{'1-NN':>7s}{'LDA':>7s}")
    for m in feats["train"][0]:
        Xtr = np.array([f[m] for f in feats["train"]])
        Xte = np.array([f[m] for f in feats["test"]])
        acc = classify(Xtr, np.array(ys["train"]), Xte, np.array(ys["test"]))
        print(f"  {m:26s}{acc[0]:7.3f}{acc[1]:7.3f}")

C1 = [("gaussian", [ell], [a]) for ell in (2.5, 3, 3.6) for a in (1, 1.5)]
experiment("C1 (Gaussian fields; classes: correlation length 2.5, 3, 3.6 x anisotropy 1, 1.5)", C1, 40)
C2 = [(kind, [2, 3, 4.5], [1, 2.5]) for kind in ("gaussian", "shot")]
experiment("C2 (Gaussian against shot noise with the same spectrum; length and anisotropy drawn at random)", C2, 60)

Output:

C1 (Gaussian fields; classes: correlation length 2.5, 3, 3.6 x anisotropy 1, 1.5): 6 classes, 40 training and 40 test images per class, random orientation
  descriptor                   1-NN    LDA
  grey histogram (16)         0.171  0.183
  GLCM directional (48)       1.000  0.921
  GLCM sorted (48)            0.996  0.992
  Gabor directional (18)      1.000  0.996
  Gabor sorted (18)           1.000  1.000
  Gabor sorted + CV (36)      0.967  1.000
  LBP riu2 8,1+16,2 (28)      0.996  1.000

C2 (Gaussian against shot noise with the same spectrum; length and anisotropy drawn at random): 2 classes, 60 training and 60 test images per class, random orientation
  descriptor                   1-NN    LDA
  grey histogram (16)         0.900  0.925
  GLCM directional (48)       0.783  0.825
  GLCM sorted (48)            0.800  0.842
  Gabor directional (18)      0.442  0.483
  Gabor sorted (18)           0.475  0.525
  Gabor sorted + CV (36)      0.900  0.933
  LBP riu2 8,1+16,2 (28)      0.900  0.950

C1: second-order structure, not grey levels. All Gaussian classes have the same standard normal grey-level distribution, so the grey-level histogram is at chance (0.171 and 0.183 against $1/6 = 0.167$). Every other descriptor separates the six classes almost perfectly (0.921–1.000), even with ratios of 1.2 between neighbouring correlation lengths: in our reading, a $128 \times 128$ image contains enough correlation lengths for the averages to be stable. The lowest value is the directional GLCM with LDA (0.921, against 0.992 for its sorted form): with random orientations, the directional features of one class vary with the angle, which, in our reading, a linear boundary absorbs less well than a nearest-neighbour rule, while sorting removes most of the variation.

C2: the same spectrum. Mean Gabor energies, directional or sorted, are at chance (0.442–0.525), as the derivation requires: whatever their parameters, the two kinds have the same expected energy in every filter. Adding the coefficient of variation of the energy, a fourth-order statistic, raises the same filter bank to 0.900 and 0.933. The grey-level histogram separates the kinds (0.900, 0.925) through the skewed distribution of the shot noise; the GLCM features reach 0.783–0.842 and the LBP histograms 0.900 and 0.950. Every planned outcome of Part C held: in C1 the grey-level histogram stayed at or below 0.30 and every other descriptor reached at least 0.90; in C2 the mean-energy descriptors stayed at or below 0.65 and every other descriptor reached at least 0.70, with both classifiers. The example shows, for classifiers, that the power spectrum, and with it any mean filter energy, is blind to a difference that other statistics see at once. It is not a test of the second-order Julesz conjecture: the two kinds differ in their grey-level histograms and in their pair distributions, which the grey-level histogram and the GLCM detect, whereas the counterexamples Portilla and Simoncelli describe have identical statistics through second order.

Common Misconceptions

• “Co-occurrence matrices are the same as the autocorrelation.” The GLCM is the full joint distribution of grey-level pairs at an offset; the autocorrelation is only its correlation. They carry the same information for Gaussian fields, as Part B shows, and different information otherwise.
• “A filter bank captures the texture.” The mean energy in each band is a function of the power spectrum alone (derived above). Banks discriminate beyond the spectrum only through nonlinear statistics of the responses, such as variances, histograms, texton assignments or the magnitude correlations of Portilla and Simoncelli’s model.
• “Texture is high-frequency detail.” Texture statistics span all scales, and the correlation lengths in the example range over a factor of four; a smooth, coarse texture is still a texture.
• “More statistics are always better.” Portilla and Simoncelli’s texture page argues that a representation is only interesting when it is minimal: with very many constraints, more than the number of input pixels, only the original image and its translations satisfy them, and synthesis reproduces the example instead of producing a new sample.
• “scikit-image’s angle $\pi/2$ looks upward.” In version 0.25.2 it pairs each pixel with the one below it; see Properties.

Where It Is Used

Material and texture recognition. Cimpoi et al. (2014) compared local descriptors on their Describable Textures Dataset (DTD), 5,640 images annotated with 47 attributes, all encoded as a bag of visual words with k-means, in mean average precision of the 47 attribute classifiers averaged over ten splits (Table 1, read from the rendered page). With their best kernel (exponential $\chi^2$), the Leung–Malik filter bank reached 39.7%, the MR8 filter bank 30.7%, raw $3 \times 3$ and $7 \times 7$ patches 30.7% and 37.1%, uniform LBP 24.8% and vector-quantized LBP 34.7%, against 53.8% for dense SIFT, a descriptor designed for objects. Cimpoi, Maji and Vedaldi (2015) then treated the convolutional layers of an ImageNet-trained network as a filter bank and pooled their outputs orderlessly with a Fisher vector (FV-CNN). In their Table 1 (read from the rendered page), on DTD, with the region equal to the whole image, accuracy was 58.6% for the Fisher vector of dense SIFT, 72.3% for FV-CNN on the 19-layer VGG-VD network against 62.9% for its fully connected layer, and 74.7% for the two combined; on the Flickr Material Dataset the corresponding values were 59.8%, 79.8%, 77.4% and 82.4%. The paper’s text gives slightly different values for some comparisons than its table (for example 66.6% for FV-CNN on DTD where the table prints 66.8% for FV-CNN on VGG-M, and 59.2% for dense SIFT on FMD where it prints 59.8%); the values here are the table’s.

Segmentation. Malik et al. (2001) segmented images with normalized cuts on affinities that compare texton histograms in local windows, combined with contour cues. TextonBoost, as the semantic segmentation article describes, classified pixels with boosted texton-based shape and texture features.

Synthesis. Portilla and Simoncelli (2000) synthesize a new sample of a texture by projecting white noise, step by step, onto the set of images that share the example’s statistics. Neural style transfer grew out of texture synthesis from an example with the Gram matrices of a convolutional network’s features as the statistics, as that article describes.

Scenes, inspection and counting. GIST describes a whole photograph by the magnitudes of Gabor-like filter responses averaged over a coarse grid (see scene recognition); industrial inspection uses statistical, structural, filter-based and model-based texture descriptors to flag regions of a regular texture such as fabric whose texture departs from the normal pattern (see industrial visual inspection); and crowd counting by regression has used the grey-level co-occurrence texture of each segment as a feature, as the object counting article reports.