# Tests of isotropy for rough textures of trended images

Frédéric Jp Richard

► To cite this version:

Frédéric Jp Richard. Tests of isotropy for rough textures of trended images. Statistica Sinica, 2016, 26, pp.1279-1304. 10.5705/ss.202014.0077. hal-01075870v2

HAL Id: hal-01075870

https://hal.science/hal-01075870v2

Submitted on 7 Sep 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tests of isotropy for rough textures of trended images.

Frédéric J.P. RICHARD

*Abstract:*

In this paper, we propose a statistical methodology to test whether the texture of an image is isotropic or not. This methodology is based on the well-known quadratic variations defined as averages of square image increments. Specific to our approach, these variations are computed in different directions using grid-preserving image rotations. We study these variations asymptotically in a framework of intrinsic random fields allowing us to take into account the presence of polynomial trends in images. We establish a convergence result linking variation and scale logarithms through an asymptotic Gaussian linear model. This model involves direction-dependent intercepts which are equal when textures are isotropic. Hence, we test the texture isotropy using Fisher tests that check the validity of the assumption of the intercept equality. These tests are validated using 6000 realizations of anisotropic fractional Brownian fields simulated on a grid of size 100 × 100. Results show that more than 70% of anisotropic cases can be detected with less than 1% of misclassified isotropic cases. Eventually, we apply our methodology to the segmentation of some biological microscopic images.

*Key words and phrases:* texture analysis, image processing, isotropy, anisotropy, intrinsic random field, statistical test, fractional Brownian field, quadratic variations, biological imaging, microscopy.

## 1 Introduction

The notion of Texture usually refers to an image aspect which is essential for processing images (see Davies (2008); Tuceryan and Jain (1998) and references therein for an introduction to Texture). From one study to another, the definition of textures varies depending mainly on analysis approaches, those approaches being based on statistical models (Besag (1986); Cross and Jain (1983); Geman and Geman (1984); Kaplan and Kuo (1995); Pentland (1984)), variational image
<figure><img src="image_2.png" /><figcaption>Figure 1: Isotropy and anisotropy of textures of images with trends. (a.1) Realization $X_a$ of a (isotropic) fractional Brownian field with constant Hurst function $\beta \equiv 0.22$, (b.1) realization $X_b$ of a anisotropic fractional Brownian field with a Hurst function ranging in (0.13, 0.76), (a.2) $X_a + P_1$, where $P_1$ is an isotropic polynomial of degree 2, (b.2) $X_b + P_1$, (a.3) $X_a + P_2$, where $P_2$ is an anisotropic polynomial of degree 2, (b.3) $X_b + P_2$.</figcaption></figure>

decompositions (Aujol et al. (2006, 2005); Meyer (2001)), patches (Efros and Leung (1999); Peyré (2009); Wei et al. (2009)), image filters (Demanet and Ying (2007); Jain and Farrokhsia (1991); Portilla and Simoncelli (2000); Unser (1995); Zhu et al. (1998)), textons (Desolneux et al. (2012); Galerne et al. (2011); Haralick (1979)), etc. In this work, we deal with rough textures that appear on irregular images (this irregularity will be given a Holderian meaning in Section 2.2). From this viewpoint, textures are associated to high frequency components of images, and not related to their possible trends and other low frequency components. This is illustrated in Figure 1 where some images having textures of interest are shown with and without trend.

In these examples, we can further distinguish two types of textures: the
isotropic and anisotropic ones. Isotropic textures (first row) are uniform in all directions whereas the anisotropic ones (second row) are oriented. Anisotropy is one of the main features for the analysis of textures. It was early recognized as an important feature for pre-attentive and attentive vision (Julesz and Bergen (1983); Rao and Lohse (1993)). It is also relevant for the characterization of material properties in Material Science (Germain et al. (2003)), the diagnostic or prognostic of diseases in Medicine (Biermé et al. (2009); Biermé and Richard (2011); Brunet-Imbault et al. (2005); Lemineur et al. (2004); Richard (2015); Richard and Biermé (2010)), the fingerprint authentication in Biometry (Jain et al. (1997); Jiang (2005)), the analysis of chemical movements through aquifers in Hydro-Geology (Benson et al. (2006)), etc.

In this paper, our goal is to construct statistical tests to decide whether the texture of an irregular image is isotropic or not. The tests are focused on the texture of the image, and do not concern its trend. In particular, images of the first row of Figure 1 should all be considered as isotropic, even though the third one has an anisotropic trend.

Richard and Biermé (2010) addressed the issue of testing the texture isotropy in the framework of anisotropic fractional Brownian fields (AFBF). Introduced by Bonami and Estrade (2003), AFBF are some Gaussian random fields that are spatial extensions on $\mathbb{R}^d$ ($d > 1$) of the fractional Brownian motion (Mandelbrot and Van Ness (1968)). Having zero-mean and stationary increments, they are characterized by a semi-variogram $v_Z$ whose harmonizable representation

$$ \forall x, y \in \mathbb{R}^d, v_Z(x) = \frac{1}{2}\mathbb{E}((Z(x + y) - Z(y))^2) = \int_{\mathbb{R}^d} |e^{i\langle w, x \rangle} - 1|^2 f(w)dw \quad (1) $$

is specified by a spectral density $f$ of the form

$$ f(w) = \tau \left( \frac{w}{|w|} \right) |w|^{-2\beta \left( \frac{w}{|w|} \right) - d}. \quad (2) $$

In this expression, the so-called topothesy function $\tau$ and Hurst function $\beta$ are two even, positive, and bounded functions defined on the unit sphere $S^{d-1} = \{w \in \mathbb{R}^d, |w| = 1\}$ of $\mathbb{R}^d$. They are two directional functions that make the model anisotropic. In particular cases when these functions are constant, the model becomes isotropic. The Hurst function is further assumed to range in an interval $(H, \overline{H}) \subset (0, 1)$ for almost points of $S^{d-1}$. An AFBF has irregular realizations
(see Section 2.2 for details). Hence, due to its properties, AFBF are well-suited to model anisotropic textures of irregular images. Some realizations of AFBF (in isotropic and anisotropic cases) are shown in Figure 1. Richard and Biermé (2010) built an isotropy test upon the statistical validation of an hypothesis of the form $\beta(s_1) = \beta(s_2) = H$ for a couple $(s_1, s_2)$ of directions. This hypothesis only gives a necessary condition for the field isotropy, and is focused on two chosen directions. Hence, the test is likely to fail detecting a texture anisotropy which checks this condition. Moreover, the test does not cover anisotropies induced by the topothesy function.

Besides, the AFBF model relies upon some stationarity assumptions which restrict their application to images without any trend. To overcome this limitation, we propose to work in a general framework of Gaussian intrinsic random fields (IRF), which is already very popular in Spatial Statistic (see Chilès and Delfiner (2012); Cressie (1993); Matheron (1973)). These fields generalize fields with stationary increments. They are characterized by a so-called generalized covariance which can be represented in an harmonizable form extending the one in Equation (1) (see Gelfand and Vilenkin (1964); Matheron (1973) and Section 2.1 for details). Using this representation, we focus on IRF having a spectral density $f$ which satisfies

$$|w| > A \Rightarrow f(w) - \tau \left( \frac{w}{|w|} \right) |w|^{-2\beta\left(\frac{w}{|w|}\right)} \le C|w|^{-2H-d-\gamma}. \quad (3)$$

for some positive constants $A, C$ and $\gamma$, and some topothesy and Hurst functions $\tau$ and $\beta$. Such a condition implies that the field texture is similar to the one of an AFBF. However, this condition only constrains high frequencies of the field. At low frequency, the spectral density only fulfills general requirements of an IRF density, allowing the field to have a polynomial trend of an arbitrary order. In this work, images are considered as realizations of Gaussian IRF satisfying (3). The texture of such an image is then called isotropic whenever

$$\tau(s) = \tau_0 \in \mathbb{R}^+, \beta(s) = H \in (0, 1) \quad (4)$$

holds almost everywhere on $S^{d-1}$. We construct statistical tests to validate the hypothesis that a given image is the realization of an isotropic IRF satisfying (4). As opposed to the one of Richard and Biermé (2010), these tests are designed to detect texture anisotropies whatever their direction and type.
The construction of our test statistics is based on quadratic variations. These variations are averages of square increments computed on an observation grid (see Section 2.1 for a definition of increments). In the literature, they are commonly used to estimate parameters characterizing a field irregularity, e.g. the Hurst index of fields related to the fractional Brownian motion (Chan and Wood (2000); Davies and Hall (1999); Kent and Wood (1997); Zhu and Stein (2002)), or the local Hurst function of the multifractional Brownian fields and their generalizations (Ayache et al. (2007); Ayache and Lévy-Véhel (2004); Benassi et al. (2000, 1998); Coeurjolly (2005); Istas and Lang (1997)). Biermé and Richard (2008); Richard and Biermé (2010) also used such variations in combination with a Radon transform for the estimation of the Hurst function of an AFBF. Due to a discretization issue of the Radon transform, their application was limited to the estimation of function values in vertical and horizontal directions. Variations we use have some original features. First, they are computed on rotated and rescaled images so as to provide both direction and scale information. Image rotations and rescalings are carefully chosen to preserve the observation grid and avoid data interpo- lations. Furthermore, the order of increments is adapted to the order of the polynomial trend so as to annihilate its effect on the estimation. We establish a convergence result linking variation and scale logarithms through an asymptotic Gaussian linear model. This model includes some direction-dependent intercepts that are equal when the texture isotropic. Hence, we propose to test the texture isotropy using some Fisher tests that verify the assumption of the equality of these intercepts.

## 2 Image and Texture Modeling

In this work, an image is considered as a realization of a random field. Let $d > 1$ be the dimension of the image (usually, $d = 2$ or $d = 3$), and $Z$ be a random field defined on $\mathbb{R}^d$, a set of random variables defined on a probability space $(\Omega, \mathcal{A}, \mathbb{P})$, and indexed in $\mathbb{R}^d$. For $N \in \mathbb{N}^*$, $k \in \mathbb{Z}^d$ and $x \in \mathbb{R}^d$, $Z(x)$ and $Z^N[k]$ denote random variables of $Z$ at positions $x$ and $k/N$, respectively. An image is defined as a realization of random variables $Z_N[k]$ at points $k$ of the grid $G_N = \llbracket 1, N \rrbracket^d$. We assume that the field $Z$ is Gaussian, i.e. any uplet $(Z(x_1), \dots, Z(x_n))$ is a
Gaussian vector. We also assume that Z is mean square continuous at any point of $\mathbb{R}^d$.

In what follows, we recall some elements about intrinsic random fields (IRF) which are required to set our framework for the modeling of random fields (interested readers may refer to Biermé (2005); Chilès and Delfiner (2012); Cressie (1993); Matheron (1973); Yaglom (1986) for more comprehensive presentations of IRF).

## 2.1 Intrinsic random fields

In the sequel, we use some shorthands for manipulating multivariate polynomials: given a multi-index $l = (l_1, \dots, l_d) \in \mathbb{N}^d$, we denote by $x^l$ the monomial $x_1^{l_1} \cdots x_d^{l_d}$ defined for $x = (x_1, \dots, x_d) \in \mathbb{R}^d$, and by $|l| = \sum_{j=1}^d l_j$ its order.

**Definition 2.1.** Let Z be a random field on $\mathbb{R}^d$ and $M \in \mathbb{N}$. An increment of order M (or M-increment) of Z is a linear combination $Z_{\lambda,x}$ of random variables from Z

$$Z_{\lambda,x} = \sum_{i=1}^{m} \lambda_i Z(x_i), \qquad (5)$$

defined with a set $\lambda = (\lambda_i)_{i=1}^m$ of scalar values and a set $x = (x_i)_{i=1}^m$ of points in $\mathbb{R}^d$ satisfying the condition

$$\sum_{i=1}^{m} \lambda_i x_i^l = 0, \forall l \in \mathbb{N}^d, |l| \le M. \qquad (6)$$

Due to (6), the linear combination defined by $\lambda$ and $x$ removes from Z any polynomial trend of order less or equal to M. Random variables $Z(x_i) - Z(x_j)$, defined with arbitrary $x_i \ne x_j$, are examples of 0-increments. In Section 3.1, we construct increments of arbitrary orders on a grid.

**Definition 2.2.** Let $M \in \mathbb{N}$. A random field Z is an intrinsic random field of order M (or M-IRF) if, for any sets $\lambda = (\lambda_i)_{i=1}^m \in \mathbb{R}^m$ and $x = (x_i)_{i=1}^m$ of points in $\mathbb{R}^d$ satisfying (6), the random field defined for all $y \in \mathbb{R}^d$ by

$$V_{\lambda,x}(y) = \sum_{i=1}^{m} \lambda_i Z(x_i + y),$$
has zero mean, and is second-order stationary,

$$ \mathbb{E}(V_{\lambda,x}(y)) = 0, \forall y \in \mathbb{R}^d, $$

$$ \mathbb{E}(V_{\lambda,x}(y)V_{\lambda,x}(z)) = K_{\lambda,x}(y-z), \forall y, z \in \mathbb{R}^d. $$

Field models previously used by Biermé et al. (2011); Chan and Wood (2000); Richard and Biermé (2010); Zhu and Stein (2002) are continuous 0-IRF. Since their increments have zero mean, their expectation is only a constant. Consequently, they cannot be used to describe fields having a polynomial trend of degree higher than 0. Such a description becomes possible with IRF of higher orders.

**Example 2.1.** Let $W$ be a continuous zero-mean 0-IRF and $(\alpha_l)_l$ a set of square integrable random variables, the random field

$$ Z(x) = \sum_{l \in \mathbb{N}^d, |l| \le M} \alpha_l x^l + W(x), $$

is a continuous $M$-IRF with a polynomial trend of degree $M$.

The correlation structure of a continuous $M$-IRF $Z$ is characterized by a generalized covariance, defined as a function $K_Z$ satisfying

$$ \mathbb{E}(Z_{\lambda,x} Z_{\mu,y}) = \sum_{i=1}^{m} \sum_{j=1}^{n} \lambda_i \mu_j K_Z(x_i - y_j) \quad (7) $$

for any pair of $M$-increments $Z_{\lambda,x}$ and $Z_{\mu,y}$ of $Z$. Generalized covariances extend ordinary covariances of stationary fields, for which (7) holds for any linear combination of $Z$. Besides, when $Z$ is a 0-IRF with a generalized covariance $K_Z$, its semi-variogram $v_Z(h) = -K_Z(h)$. A generalized covariance plays the same characterization role for a $M$-IRF as the semi-variogram for a 0-IRF.

For any continuous $M$-IRF $Z$, there exists a continuous generalized covariance $K_Z$, which is unique up to an even polynomial of order $2M$ (Gelfand and Vilenkin (1964); Matheron (1973)). This function is symmetric and $M$-conditionally positive definite,

$$ \sum_{i,j} \lambda_i \lambda_j K_Z(x_i - x_j) \ge 0, $$
for any sets $(\lambda_i)_i$ and $(x_i)_i$ satisfying (6). Reciprocally, any continuous, symmetric and $M$-conditionally positive definite is a generalized covariance of a continuous $M$-IRF.

Besides, $M$-conditionally positive definite functions can be characterized using a spectral theorem established by Gelfand and Vilenkin (1964) for generalized random functions, and developed specifically by Matheron (1973) for ordinary random fields. We recall a corollary of this theorem which gives a sufficient condition for a function to be a generalized covariance.

**Theorem 2.3.** A function $K$ on $\mathbb{R}^d$ is a generalized covariance of an $M$-IRF if it is of the form

$$K(h) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} (\cos(\langle w, h \rangle) - \mathbf{1}_B(w) P_M(\langle w, h \rangle)) f(w) dw + Q(h), \quad (8)$$

where $P_M(t) = 1 - \frac{t^2}{2} + \cdots + \frac{(-1)^M}{(2M)!} t^{2M}$, $\mathbf{1}_B(w)$ is the indicator function of an arbitrary neighborhood of 0, $Q$ an arbitrary even polynomial of degree $\le 2M$, and $f$ is an even and positive function satisfying the integrability conditions

$$\forall A > 0, \quad \int_{|w|<A} |w|^{2M+2} f(w) dw < \infty \quad \text{and} \quad \int_{|w|>A} f(w) dw < \infty. \quad (9)$$

The function $f$ is called the spectral density of the field. Conditions (9) concern the spectral density at low and high frequencies, respectively. In particular, the first condition becomes weaker as $M$ increases. This enables inclusion of densities having a larger power at low frequencies, and to account for fields having polynomial trends of larger order.

**Example 2.2.** The anisotropic fractional Brownian fields, whose spectral density has the form (2), are examples of 0-IRF. Such fields can be easily extended to arbitrary IRF by letting

$$\overline{H} = \operatorname{ess\,sup}\{\beta(s), s \in S^{d-1}, \tau(s) \neq 0\} \quad (10)$$

vary in $(0, +\infty)$ instead of $(0, 1)$. The obtained IRF is then of order $M \ge \overline{H} - 1$.

## 2.2 Hölder irregularity

Images we consider in this work are irregular. We define this irregularity in a Hölder sense (Adler (1981); Biermé (2005)).
**Definition 2.4.** A field Z satisfies a uniform stochastic Hölder condition of order $\alpha \in (0, 1)$ if, for any compact set $C \subset \mathbb{R}^d$, there exists an almost surely finite, positive random variable $A_C$ such that the Hölder condition

$$|Z(x) - Z(y)| \le A_C |x - y|^\alpha. \qquad (11)$$

holds for any $x, y \in C$, with probability one.

If there exists $H \in (0, 1)$ for which (11) holds for any $\alpha < H$ but not for $\alpha > H$, then we say that $Z$ admits $H$ as critical Hölder exponent, or that $Z$ is $H$-Hölder.

The critical Hölder exponent $H$ of a field $Z$ quantifies its irregularity. It is close to 0 for the most irregular fields and to 1 for the least ones. Visually, field textures look rougher and rougher as $H$ decreases to 0.

The Hölder irregularity of a Gaussian M-IRF with spectral density can be characterized from the asymptotic behavior of its spectral density at high frequency (as $|w| \to +\infty$). To be more specific, we state a proposition borrowed from Biermé (2005, Propositions 2.1.6 and 2.1.7).

**Proposition 2.5.** Let $Z$ be a continuous Gaussian IRF of an arbitrary order with a spectral density $f$, and $H \in (0, 1)$. $Z$ is $H$-Hölder if and only if, for any $0 < \alpha < H < \beta < 1$, there exist positive constants $A, B_1, B_2$, and a positive measure subset $E$ of the unit sphere $S^{d-1}$ of $\mathbb{R}^d$ such that for almost all $w \in \mathbb{R}^d$

$$|w| \ge A \Rightarrow f(w) \le B_2|w|^{-2\alpha-d}, \qquad (12)$$

$$|w| \ge A \text{ and } \frac{w}{|w|} \in E \Rightarrow f(w) \ge B_1|w|^{-2\beta-d}. \qquad (13)$$

Condition (12) (resp. (13)) ensures that the critical Hölder exponent of a Gaussian IRF is above (resp. below) $H$.

**Example 2.3.** In the particular case of an anisotropic fractional Brownian field of density (2), Conditions (12) and (13) are fulfilled for

$$H = \operatorname{ess\,inf}\{\beta(s), s \in S^{d-1}, \tau(s) \neq 0\}. \qquad (14)$$

It is the same for an IRF whose spectral density satisfies (3). Hence, both fields are $H$-Hölder.
# 3 Test construction

## 3.1 Quadratic variations

First, we describe increments (see Definition 2.1) that compose quadratic variations. These increments are computed by application to the discrete field $Z^N$ of a linear filter

$$ \forall m \in \mathbb{Z}^d, V^N[m] = \sum_{k \in [0, L]^d} v[k] Z^N[m-k], \quad (15) $$

determined by a convolution kernel $v$ of a finite support $[0, L]^d$ with $L \in \mathbb{N}^d$. Let $Q_v$ be the characteristic polynomial associated to $v$

$$ \forall z \in \mathbb{R}^d, Q_v(z) = \sum_{k \in [0, L]^d} v[k] z^k. \quad (16) $$

We now state a necessary and sufficient condition on partial derivatives of $Q$ for the random variable $V^N[m]$ to be a $K$-increment of $Z$ for any $m \in \mathbb{Z}^d$ (see Section S2 of Supplementary Material for a proof).

**Proposition 3.1.** Let $K \in \mathbb{N}$. For all $m \in \mathbb{Z}^d$, random variables $V^N[m]$ are $K$-increments of $Z$ if and only if

$$ \forall l \in [0, K]^d, |l| \le K, \frac{\partial^{|l|} Q_v}{\partial z_1^{l_1} \dots \partial z_d^{l_d}}(1, \dots, 1) = 0. \quad (17) $$

**Example 3.1.** Let $L \in \mathbb{N}^d$, $|L| > 0$, and $v$ be the kernel associated to the characteristic polynomial $Q_v(z) = (z_1 - 1)^{L_1} \dots (z_d - 1)^{L_d}$. From Proposition 3.1, it follows that $v$ induces increments of order $K = |L| - 1$. Its terms are

$$ v[l] = (-1)^{|l|} \binom{L_1}{l_1} \dots \binom{L_d}{l_d}, \quad \text{if } l \in [0, L]^d, v[l] = 0 \text{ otherwise,} \quad (18) $$

where $\binom{n}{k}$ stands for the binomial coefficient. For $0 < |L| \le 2$ and $d = 2$, such kernels yield the main increments used by Biermé et al. (2011); Chan and Wood (2000); Zhu and Stein (2002) for the estimation of Hölder irregularity.

Our statistical test involves computing increments at different scales and orientations. To do so, the observation field $Z^N$ is transformed into discrete
fields $Z_T^N$ by applying combinations $T$ of rescaling and rotation that map $\mathbb{Z}^d$ into itself:

$$ \forall k \in \mathbb{Z}^d, Z_T^N[k] = Z\left(\frac{Tk}{N}\right) = Z^N[Tk]. \qquad (19) $$

When $d=2$, such transforms have the form

$$ T_u = \begin{pmatrix} u_1 & -u_2 \\ u_2 & u_1 \end{pmatrix} = |u| \begin{pmatrix} \cos(\arg(u)) & -\sin(\arg(u)) \\ \sin(\arg(u)) & \cos(\arg(u)) \end{pmatrix}, \quad \text{for } u \in \mathbb{Z}^2 \setminus \{(0,0)\}, $$

corresponding to a rescaling of factor $|u|$ and a rotation of angle $\arg(u)$. When $d=3$, examples of transforms are

$$ T_u = \begin{pmatrix} u_1 & -u_2 & 0 \\ u_2 & u_1 & 0 \\ 0 & 0 & |u| \end{pmatrix}, \begin{pmatrix} u_1 & 0 & -u_2 \\ 0 & |u| & 0 \\ u_2 & 0 & u_1 \end{pmatrix}, \begin{pmatrix} |u| & 0 & 0 \\ 0 & u_1 & -u_2 \\ 0 & u_2 & u_1 \end{pmatrix}. $$

Transforms can be identified by a single vector $u$ in $\mathbb{R}^2$ or $\mathbb{R}^3$, respectively. The application of a filter $v$ to a transform field $Z_T^N$ leads to random variables

$$ \forall m \in \mathbb{Z}^d, V_T^N[m] = \sum_k v[k] Z\left(\frac{m-Tk}{N}\right) = \sum_k v[k] Z^N[m-Tk]. \qquad (20) $$

**Proposition 3.2.** Random variables $V_T^N[m]$ are K-increments of Z if and only if the characteristic polynomial of $v$ satisfies (17).

A proof is given in Section S2 of Supplementary Material.

Next, we select a set of transforms. We use a multi-index $a = (i, j)$ to classify transforms according to their rotation angle: the transform $T_a$ is identified by a vector $u_a$. It corresponds to a $j$th transform among those that have a same rotation angle $\arg u_{i1}$. The rescaling factor of this transform is $|u_a|$. The index $a$ varies in the set

$$ \mathcal{F} = \{i \in [1, n_b], j \in [1, P_i]\}, $$

where $n_b$ is the number of different rotation angles and $P_i \in \mathbb{N}^*$ depends on $i$. We denote by $n_f$ the cardinality of $\mathcal{F}$.

Applying a chosen filter $v$ to transformed fields $T_a$, we obtain a set of increments denoted by $V_a^N[m]$, for $m \in \mathbb{Z}^d$. We define the vector-valued random field $V^N = ((V_a^N[m])_{a \in \mathcal{F}}, m \in \mathbb{Z}^d)$, and state some of its properties under the assumption that Z is an IRF (see Section S2 of Supplementary Materials for a proof).
**Proposition 3.3.** For some $M \in \mathbb{N}$, let $Z$ be a continuous M-IRF with a spectral density $f$. Assume that the vector-valued field $V^N$ is constructed using increments of order $K \ge M$. Then, $V^N$ is stationary with zero mean. Moreover, for any $a, b \in \mathcal{F}$ and $m, n \in \mathbb{Z}^d$,

$$ \mathbb{E}(V_a^N[m]V_b^N[n]) = \frac{1}{(2\pi)^d} \int_{[0,2\pi]^d} f_{a,b}^N(w)e^{i\langle m-n,w \rangle}dw, \quad (21) $$

where $f_{a,b}^N$, the multivariate spectral density of $V^N$, is a function of $L^1([0, 2\pi]^d)$ given by

$$ f_{a,b}^N(w) = N^d \hat{v}(T_a'w) \bar{v}(T_b'w) \sum_{k \in \mathbb{Z}^d} f(N(w + 2k\pi)), \quad (22) $$

with $T'$ the transpose of $T$, $\hat{v}(w) = \sum_{k \in \mathbb{Z}^d} v[k]e^{-i\langle k, w \rangle}$ the discrete Fourier transform of $v$, and $\bar{v}$ its conjugate.

We now consider the set $\mathcal{E}_N$ of cardinality $N_e$ including all points $m$ of $\mathbb{Z}^d$ for which all increments $V_a^N[m]$ can be computed using exclusively observations on the grid $\mathcal{G}_N$. We define the quadratic variations on $\mathcal{E}_N$ as

$$ \forall a \in \mathcal{F}, W_a^N = \frac{1}{N_e} \sum_{m \in \mathcal{E}_N} (V_a^N[m])^2. \quad (23) $$

These variations account for those of Biermé et al. (2011); Chan and Wood (2000), which can be obtained for $d = 2$ with $u \propto (1,0)$ (i.e. without any field rotation). They also contains variations defined by Zhu and Stein (2002) using some field rotations. However, they are more generic as they include variations obtained with arbitrary grid-preserving rotations and K-order filters.

For $a \in \mathcal{F}$, we also define log-variations $Y_a^N$ and log-scales $x_a$ as

$$ Y_a^N = \log(W_a^N) \quad \text{and} \quad x_a = \log(|u_a|^2), \text{ respectively.} \quad (24) $$

## 3.2 Gaussian linear model

The construction of the statistical test is based on asymptotic properties of the random vector $Y^N = (Y_a^N)_{a \in \mathcal{F}}$.

**Theorem 3.4.** Assume that $Z$ is a continuous Gaussian M-IRF with a spectral density $f$ satisfying (3). Let

$$ H = \operatorname{ess\,inf}\{\beta(s), s \in S^{d-1}, \tau(s) \neq 0\}. \quad (25) $$
Assume that $H \in (0,1)$, and the set $E_0 = \{s \in S^{d-1}, \beta(s) = H, \tau(s) \neq 0\}$ is of positive measure on $S^{d-1}$. If a log-variation vector $Y^N$ is constructed using increments of orders $K \ge M+1$, and $K \ge M/2 + d/4$ if $d > 4$, then

$$N^{\frac{d}{2}}(Y^N - \zeta^N) \xrightarrow[N \to +\infty]{d} \mathcal{N}(0, \Sigma), \quad (26)$$

for some vectors $\zeta^N \in \mathbb{R}^{n_f}$ and $a n_f \times n_f$-covariance matrix $\Sigma$. Here, for $a, b \in \mathcal{F}$, terms of the matrix $\Sigma$ are given by

$$\Sigma_{a,b} = \frac{2(2\pi)^d \int_{[0,2\pi]^d} |f_{a,b}(w)|^2 dw}{\int_{[0,2\pi]^d} f_{a,a}(w) dw \int_{[0,2\pi]^d} f_{b,b}(w) dw}, \quad (27)$$

where $f_{a,b}$ is a function defined for almost $w \in [0, 2\pi]^d$ by

$$f_{a,b}(w) = \hat{v}(T'_{a}w)\bar{v}(T'_{b}w) \sum_{k \in \mathbb{Z}^d} \delta\left(\frac{w + 2k\pi}{|w + 2k\pi|}\right) |w + 2k\pi|^{-2H-d} \quad (28)$$

with $\delta$ defined, for almost $s \in S^{d-1}$, by $\delta(s) = \lim_{\rho \to +\infty} f(s\rho)\rho^{2H+d}$. Terms $\zeta^N$ are of the form

$$\forall a = (i, j) \in \mathcal{F}, \zeta_{ij}^N = x_{ij}H + \beta_i^N, \quad (29)$$

with

$$\beta_i^N = \log(C_H(\arg(u_{i1}), v)) - 2H \log(N) - d \log(2\pi), \quad (30)$$

where

$$C_H(\arg(u_a), v) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \left| \hat{v} \left( \frac{T'_a}{|u_a|} w \right) \right|^2 \delta \left( \frac{w}{|w|} \right) |w|^{-2H-d} dw. \quad (31)$$

Moreover, when $Z$ has an isotropic texture, $\beta_i^N = \tilde{\beta}^N$ for all $i, j$.

This theorem is proved in Section S1 of Supplementary Material. It brings out a linear asymptotic dependency between variations and parameters related to field irregularity and anisotropy. This dependency can be formalized in terms of a Gaussian linear model

$$\forall (i, j) \in \mathcal{F}, Y_{ij}^N = x_{ij}H + \beta_i^N + \epsilon_{ij}^N, \quad (32)$$

where $\epsilon_{ij}^N$ are correlated Gaussian variables. It can be interpreted as a generalized model of analysis of covariance, in which the $x_{ij}$ are regression variables, and the
rotation angle (index j) used to form $Y_{ij}^N$ is a qualitative factor. As mentioned in Theorem 3.4, parameters $\beta_i$ vary according to the rotation index, and are constant when the texture is isotropic. Hence, for testing the texture isotropy, we propose to check statistically the null hypothesis

$$ \mathcal{H}_0: \forall i \in [1, n_b], \beta_i = \tilde{\beta}, \quad (\text{isotropy}) \qquad (33) $$

against its alternative hypothesis $\mathcal{H}_1: \exists i, j, \beta_i \neq \beta_j$ (anisotropy).

### 3.3 Parameter estimation

Model (32) can also be formulated in a matricial form. Let $Y^N$ and $\epsilon^N$ be random vectors of size $n_f$ formed by terms $Y_a^N$ and $\epsilon_a^N$, respectively. Define a parameter vector $\theta$ of size $n_b+1$. Set a design matrix $X$, of size $n_f \times (n_b+1)$, having terms $X_{ij,0} = x_{ij}$ on the $(i, j)^{th}$ row and the first column and, $X_{ij,m} = \delta_{i=m}$ on $(i, j)^{th}$ row of the $m^{th}$ column. Then, (32) is equivalent to

$$ Y^N = X\theta + \epsilon^N, \epsilon^N \sim \mathcal{N}(0, \Sigma). \qquad (34) $$

Using a generalized least square criterion, parameters $\theta$ of this model can be estimated by

$$ \tilde{\theta}^N = P(\Sigma)Y^N \quad \text{with} \quad P(\Sigma) = (X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}. \qquad (35) $$

In particular, the field irregularity $H$ is estimated by

$$ \tilde{H}^N = \hat{\theta}_0^N = P_0(\Sigma)Y^N, \qquad (36) $$

$P_0(\Sigma)$ being the first row of the matrix $P(\Sigma)$. Similarly, for $i \in [1, n_b]$, parameters $\beta_i^N$ are estimated by $\tilde{\theta}_i^N$.

### 3.4 Fisher test

To test hypothesis $\mathcal{H}_0$ against $\mathcal{H}_1$, we construct a Fisher test. Under $\mathcal{H}_0$, (32) reduces to the model

$$ Y^N = \tilde{X}\alpha + \epsilon^N, \epsilon^N \sim \mathcal{N}(0, \Sigma), \qquad (37) $$

where $\alpha$ is a parameter vector of size 2, and $\tilde{X} = (x_a)_{a \in \mathcal{F}}$ is the design matrix of size $n_f \times 2$ having terms $\tilde{X}_{ij,0} = x_{ij}$ on the $(i, j)^{th}$ row and the first column
and, 1 on rows of the second column. The GLS estimator of parameters of this model is

$$ \hat{\alpha}^N = \tilde{P}(\Sigma)Y^N \quad \text{with} \quad \tilde{P}(\Sigma) = (\tilde{X}'\Sigma^{-1}\tilde{X})^{-1}\tilde{X}'\Sigma^{-1}, \qquad (38) $$

Then, we define the test statistic

$$ F^N(\Sigma) = \frac{n_f - n_b - 1}{n_b - 1} \frac{|\Sigma^{-1/2}(X\hat{\beta}^N - \tilde{X}\hat{\alpha}^N)|^2}{|\Sigma^{-1/2}(Y - X\hat{\beta}^N)|^2}. \qquad (39) $$

Under $\mathcal{H}_0$, $F^N(\Sigma)$ has a Fisher distribution with $(n_b - 1, n_f - n_b - 1)$ as degrees of freedom. Let $\alpha \in (0, 1)$, and set $s_\alpha$ to be the quantile of order $1 - \alpha$ of a $(n_b - 1, n_f - n_b - 1)$-Fisher distribution. Then, the rejection interval $R_\alpha = \{F^N > s_\alpha\}$ defines an isotropy test with a confidence level $\alpha$.

The expression of the covariance matrix $\Sigma$ depends on the spectral density of $Z$, which is unknown. Hence, to implement the test, we compute $F^N(\tilde{\Sigma}^N)$ instead of $F^N(\tilde{\Sigma})$ using an estimate $\tilde{\Sigma}^N$ of the covariance matrix of the log-variations (see Section S3 of Supplementary Material for more details).

# 4 Numerical results

We conducted some experiments on AFBF (see (2)) simulated using the turning-band method developed by Biermé et al. (2015). This method was applied with approximately 500 bands to simulate field realizations on a grid of size $100 \times 100$. We simulated fields by triplets consisting of an isotropic fractional Brownian field (an AFBF with constant Hurst and topothesy functions), an AFBF with a non-constant Hurst function and a constant topothesy (a Hurst-induced anisotropy), and an AFBF with a constant Hurst function and a non-constant topothesy function (a topothesy-induced anisotropy). For anisotropic fields, the Hurst and topothesy functions were specified using a function $g_{a_1,a_2,\delta_1,\delta_2,\varphi}$ depending on a few parameters $a_1 > 0, a_2 > 0, 0 < \delta_1 \le \pi/2, 0 < \delta_2 < \pi/2 - \delta_1, \varphi \in (0, \pi)$. This function was defined as the $\pi$-periodic and $C^2$ function equal to $a_1$ on a set of intervals $\cup_{k \in \mathbb{Z}}(\varphi - \delta_1 + k\pi, \varphi + \delta_1 + k\pi), a_2$ in another disjoint set $\cup_{k \in \mathbb{Z}}(\varphi + \pi/2 - \delta_2 + k\pi, \varphi + \pi/2 + \delta_2 + k\pi)$, and polynomial of order 2 outside those two sets. Such a function leads to a field whose the anisotropy is oriented in a privileged direction $\varphi$ (see Examples in Figure 2). For each field triplet, we set parameters $\delta_1, \delta_2, \varphi$ by sampling from uniform distributions on $(0, \pi/2)$,
<figure><img src="image_3.png" /></figure>

$\beta \equiv 0.41, \tau \equiv 1,$
$p = 0.95.$

<figure><img src="image_4.png" /></figure>

$\beta \equiv 0.12, \tau \equiv 1,$
$p = 0.5.$

<figure><img src="image_5.png" /></figure>

$\beta \equiv 0.6, \tau \equiv 1,$
$p = 0.1.$

<figure><img src="image_6.png" /></figure>

$\beta(s) \in [0.18, 0.21],$
$\tau \equiv 1, p = 10^{-2}.$

<figure><img src="image_7.png" /></figure>

$\beta(s) \in [0.12; 0.58],$
$\tau \equiv 1, p = 10^{-8}.$

<figure><img src="image_8.png" /></figure>

$\beta(s) \in [0.33; 0.55],$
$\tau \equiv 1, p = 10^{-10}.$

<figure><img src="image_9.png" /></figure>

$\tau(s) \in [1.24, 1.52],$
$\beta \equiv 0.83, p = 10^{-2}.$

<figure><img src="image_10.png" /></figure>

$\tau(s) \in [0.3, 1.51],$
$\beta \equiv 0.32, p = 10^{-8}.$

<figure><img src="image_11.png" /></figure>

$\tau(s) \in [0.22, 0.87],$
$\beta \equiv 0.22, p = 10^{-10}.$

Figure 2: Examples of simulated fields. On the first, second, and third rows, fields are (isotropic) fractional Brownian fields, AFBF with a constant topothesy function, and AFBF with a constant Hurst function, respectively. P-values of the isotropy test (filter (0,2)) applied to each realizations are indicated by *p* below each image.
(0, π/2 − δ₁), and (0, π), respectively. We also sampled parameters H₁ ∼ U(0, 1), H₂ ∼ U(H₁, 1), C₂ ∼ |N(0, 1)|, and C₁ ∼ C₂ + |N(0, 1)|. Fields of a same triplet were then defined as an AFBF with β ≡ H₁, τ ≡ 1, an AFBF with τ ≡ 1 and β = gH₁,H₂,δ₁,δ₂,ϕ, and an AFBF with τ = gC₁,C₂,δ₁,δ₂,ϕ and β ≡ H₁. In this way, fields of a same triplet had the same order H₁ of Hölder irregularity, and anisotropic fields had a same privileged direction ϕ. We simulated 2000 triplets (6000 fields), covering largely the range of the irregularity and orientation pair (H₁, ϕ).

For testing isotropy, we used quadratic variations obtained with different increment filters of Example 3.1. We took two filters of order 1 specified ((L₁, L₂) in {(0,2), (1,1)}), two of order 2 ((L₁, L₂) in {(0,3), (1,2)}), three of order 3 ((L₁, L₂) ∈ {(0,4), (1,3), (2,2)}). According to Theorem 3.4, filters of order 1, 2 and 3 are adapted to 0-IRF, 1-IRF, and 2-IRF, respectively. We did not consider filters (L₁, L₂) with L₁ > L₂ since they are symmetric to filters (L₂, L₁) and lead to the same quadratic variations. For filters with L₁ ≠ L₂, we computed quadratic variations on images transformed by composition of rotations of angles 0 (i.e. for u ∝ (1,0)), π/4 (u ∝ (1,1)), π/2 (u ∝ (0,1)), or 3π/4 (u ∝ (1,−1)), and rescalings of factors |u|² with 1 ≤ |u| ≤ 6. Hence, we had six quadratic variations in vertical or horizontal directions (u ∝ (1,0) or u ∝ (0,1)), and four in each diagonal direction (u ∝ (1,1) or u ∝ (1,−1)). For the filter (1,1), we only used rotations of angles 0 and π/4. Indeed, due to the symmetry of filters for which L₁ = L₂, quadratic variations are the same up to a rotation of angle π/2. For each filter, we computed the Fisher statistic FN defined in (39). Given a threshold sα, a realization could then be considered as isotropic whenever FN ≤ sα, and anisotropic in the opposite case. Following Theorem 3.4, under the null assumption H₀, Fisher statistics are theoretically distributed according to a Fisher law of degree (3,15) if (L₁ ≠ L₂) and (1,7) if L₁ = L₂. Hence, the threshold sα could be set to ensure a theoretical level of significance P>H₀(FN > sα) ≤ α, for α ∈ (0, 1).

To assess the performance of a test, we estimated the probability distribution P_H1(F^N \le s_\alpha) of errors of type II. For that, we computed the ratio between the number of misclassified anisotropic cases and the total number of anisotropic case, for different significance levels α. This was done for the two anisotropy
<figure><img src="image_12.png" /><figcaption>Filter (0,2)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Filter (1,1)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Filter (0,3)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Filter (1,2)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Filter (0,4)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Filter (1,3)</figcaption></figure>

<caption>Figure 3: Evaluation of the isotropy tests on 6000 simulated fields.</caption>
types separately and jointly. We also computed the empirical distribution of errors of type I ($P_{H_0}(F^N > s_\alpha)$) to check its adequacy to the theoretical one. Evaluation results are shown in Figure 3.

We observe that empirical distributions of type I errors (black lines) are slightly above the significance level. This means that values of the Fisher statistics tend to be higher than expected on isotropic realizations. However, as can be seen in Table 1, this inflation is below 2.3% for the most commonly employed levels of confidence. For each filter, we tested the adequacy of the empirical distribution of $F^N$ (under $H_0$) to its theoretical one using a one-way Kolmogorov-Smirnov test. Adequacy assumptions were validated with p-values 0.03, 0.08, 0.3, 0.15, 0.22, 0.28, and 0.48 for filters (0, 2), (1, 1), (0, 3), (1, 2), (0, 4), (1, 3) and (2, 2), respectively.

<table><caption>Table 1: Empirical errors of type I (in %) obtained for common levels of significance.</caption><thead><tr><th>Level of sig. (%)</th><th>1</th><th>2.5</th><th>5</th><th>10</th></tr></thead><tbody><tr><td>Filter (0,2)</td><td>1</td><td>2.8</td><td>5.95</td><td>12.3</td></tr><tr><td>Filter (1,1)</td><td>0.7</td><td>2.2</td><td>4.65</td><td>9.75</td></tr><tr><td>Filter (0,3)</td><td>0.7</td><td>3.1</td><td>6.1</td><td>12.1</td></tr><tr><td>Filter (1,2)</td><td>1.1</td><td>3.4</td><td>5.9</td><td>11.7</td></tr><tr><td>Filter (0,4)</td><td>1.15</td><td>2.65</td><td>5.45</td><td>11.25</td></tr><tr><td>Filter (1,3)</td><td>0.7</td><td>2.35</td><td>5.45</td><td>11.25</td></tr></tbody></table>

Filters had different performances in terms of type II errors (red lines). The symmetric filter (1, 1) had the worst results. This filter mixes increments in two orthogonal directions, which leads to a loss of directional information. The filter (2, 2) had comparable poor performances, which are not presented in the figure. The filter (0, 3) (resp. (0, 4)) outperformed the filter (1, 2) (resp. (1, 3)) suggesting that the most directional filters of the form $(0, L_2)$ are preferable to the other ones. Choosing those filters led to fairly good results. For a significance level of 0.01, we detected 78.3%, 76.6%, and 73.3% percent of anisotropic cases for (0, 2), (0, 3) and (0, 4), respectively. The performance of these filters decreased as the filter order increased. This was probably due to the fact that the number of points ($n_e$) used for computing quadratic variations decreases as the order increases. Nevertheless, the detection score obtained with the filter (0, 4) of
higher order was still close to the filter (0, 2). Hence, we conclude that anisotropic cases can be well-detected using filters of the form (0, $L_2$) even when a polynomial trend is present in images.

We notice that percents of detected cases were significantly lower for cases with a topothesy anisotropy (red dot lines) than with the topothesy one (red dash lines). This performance difference will be explained below.

We further investigated effects of simulation and analysis parameters on the isotropic test. For that, we classified simulations of anisotropic fields into isotropic and anisotropic by applying the test of filter (0, 2) with a level of significance $\alpha = 0.05$. We then computed numbers and percentages of detected anisotropic cases by range of parameter values.

<table><caption>Table 2: Effect of the anisotropy orientation on the detection.</caption><thead><tr><th rowspan="3">Range of φ</th><th rowspan="3">Total case number</th><th colspan="4">Analysis in four directions</th><th colspan="2">Analysis in two directions</th></tr><tr><th colspan="2">Detected cases</th><th colspan="2">Detected cases</th><th rowspan="2">H. anisot.</th><th rowspan="2">T. anisot.</th></tr><tr><th>H. anisot.</th><th>T. anisot.</th><th>H. anisot.</th><th>T. anisot.</th></tr></thead><tbody><tr><td>[-1.57,-1.37]</td><td>137</td><td>125 (91%)</td><td>114 (83%)</td><td>131 (96%)</td><td>122 (89%)</td><td></td><td></td></tr><tr><td>[-1.37,-0.98]</td><td>261</td><td>231 (89%)</td><td>204 (78%)</td><td>227 (87%)</td><td>196 (75%)</td><td></td><td></td></tr><tr><td>[-0.98,-0.59]</td><td>233</td><td>211 (91%)</td><td>184 (79%)</td><td>140 (60%)</td><td>88 (38%)</td><td></td><td></td></tr><tr><td>[-0.59,-0.20]</td><td>240</td><td>214 (89%)</td><td>177 (74%)</td><td>214 (89%)</td><td>167 (70%)</td><td></td><td></td></tr><tr><td>[-0.20, 0.20]</td><td>246</td><td>220 (89%)</td><td>194 (79%)</td><td>225 (91%)</td><td>202 (82%)</td><td></td><td></td></tr><tr><td>[0.20, 0.59]</td><td>222</td><td>205 (92%)</td><td>175 (79%)</td><td>202 (91%)</td><td>164 (74%)</td><td></td><td></td></tr><tr><td>[0.59, 0.98]</td><td>289</td><td>255 (88%)</td><td>228 (79%)</td><td>165 (57%)</td><td>113 (39%)</td><td></td><td></td></tr><tr><td>[0.98, 1.37]</td><td>238</td><td>213 (89%)</td><td>195 (82%)</td><td>213 (89%)</td><td>190 (80%)</td><td></td><td></td></tr><tr><td>[1.37, 1.57]</td><td>134</td><td>120 (90%)</td><td>108 (81%)</td><td>121 (90%)</td><td>116 (87%)</td><td></td><td></td></tr></tbody></table>

We first focus on the effect of the anisotropy direction $\phi$. In Table 2, we report results obtained with two isotropy tests. The first one corresponds to the test described initially ("analysis in four directions"), whereas the second one ("analysis in two directions") is based on variations computed in the horizontal and vertical directions only. Concerning the first test, we observe that, for both anisotropy types, the number of detected cases were quite uniform with respect to the direction. This suggests that the detection is independent of this factor. This independence was confirmed by a $\chi^2$-test with a large p-value of 0.56. This clearly differed from the second test. Without the diagonal directions, the test
performances were much lower when the anisotropy orientation $\varphi$ was away from the vertical and horizontal directions (especially when $\varphi \in (-0.98, -0.59)$ or $\varphi \in (0.59, 0.98)$).

<table><caption>Table 3: Evaluation of isotropy tests of Richard and Biermé (2010).</caption><thead><tr><th rowspan="2">Range of $\varphi$</th><th rowspan="2">Total case number</th><th colspan="2">Test 1 of Richard and Biermé (2010)</th></tr><tr><th colspan="2">Detected cases</th></tr><tr><th></th><th></th><th>H. anisot.</th><th>T. anisot.</th></tr></thead><tbody><tr><td>[-1.57,-1.37]</td><td>137</td><td>21 (15%)</td><td>8 (6%)</td></tr><tr><td>[-1.37,-0.98]</td><td>261</td><td>37 (14%)</td><td>8 (3%)</td></tr><tr><td>[-0.98,-0.59]</td><td>233</td><td>41 (18%)</td><td>12 (5%)</td></tr><tr><td>[-0.59,-0.20]</td><td>240</td><td>41 (17%)</td><td>8 (3%)</td></tr><tr><td>[-0.20, 0.20]</td><td>246</td><td>47 (19%)</td><td>13 (5%)</td></tr><tr><td>[0.20, 0.59]</td><td>222</td><td>52 (23%)</td><td>14 (6%)</td></tr><tr><td>[0.59, 0.98]</td><td>289</td><td>50 (17%)</td><td>16 (6%)</td></tr><tr><td>[0.98, 1.37]</td><td>238</td><td>33 (14%)</td><td>11 (5%)</td></tr><tr><td>[1.37, 1.57]</td><td>134</td><td>19 (14%)</td><td>5 (4%)</td></tr></tbody></table><table><thead><tr><th colspan="2">Test 2 of Richard and Biermé (2010)</th></tr><tr><th colspan="2">Detected cases</th></tr><tr><th>H. anisot.</th><th>T. anisot.</th></tr></thead><tbody><tr><td>30 (22%)</td><td>9 (7%)</td></tr><tr><td>45 (17%)</td><td>12 (5%)</td></tr><tr><td>20 (9%)</td><td>4 (2%)</td></tr><tr><td>29 (12%)</td><td>21 (9%)</td></tr><tr><td>37 (15%)</td><td>22 (9%)</td></tr><tr><td>34 (15%)</td><td>21 (9%)</td></tr><tr><td>17 (6%)</td><td>18 (6%)</td></tr><tr><td>40 (17%)</td><td>15 (6%)</td></tr><tr><td>31 (23%)</td><td>4 (3%)</td></tr></tbody></table>

On the same data, we also evaluated the two isotropy tests proposed by Richard and Biermé (2010) (see Table 3). Both of these tests had much lower detection scores than the new ones. Their performances were particularly low when the anisotropy was induced by the topothesy or oriented away from the vertical and horizontal directions.

In Table 4, we deal with the order of Hölder irregularity $H_1$. We notice that
<table><caption>Table 4: Effect of the irregularity on the anisotropy detection.</caption><thead><tr><th rowspan="2">Range of H<sub>1</sub></th><th colspan="2">Detected cases</th></tr><tr><th>H. anisotropy</th><th>T. anisotropy</th></tr></thead><tbody><tr><td>[0, 0.1]</td><td>175/206 (85%)</td><td>95/206 (46%)</td></tr><tr><td>[0.1, 0.2]</td><td>197/203 (97%)</td><td>152/203 (75%)</td></tr><tr><td>[0.2, 0.3]</td><td>194/202 (96%)</td><td>160/202 (79%)</td></tr><tr><td>[0.3, 0.4]</td><td>186/194 (96%)</td><td>152/194 (78%)</td></tr><tr><td>[0.4, 0.5]</td><td>207/212 (98%)</td><td>176/212 (83%)</td></tr><tr><td>[0.5, 0.6]</td><td>186/203 (92%)</td><td>179/203 (88%)</td></tr><tr><td>[0.6, 0.7]</td><td>178/188 (95%)</td><td>161/188 (86%)</td></tr><tr><td>[0.7, 0.8]</td><td>164/181 (91%)</td><td>157/181 (87%)</td></tr><tr><td>[0.8, 0.9]</td><td>176/197 (89%)</td><td>163/197 (83%)</td></tr><tr><td>[0.9, 1]</td><td>131/214 (61%)</td><td>184/214 (86%)</td></tr></tbody></table>

the detection was worse for fields with a very low regularity ($0 < H_1 \le 0.2$) than for the more regular one. But, for a regularity above 0.2, detection performances seem quite uniform with respect to $H_1$. We applied a $\chi^2$-test on all anisotropic cases for which $H_1 > 0.2$. This test confirmed the independence of the irregularity and the detection with a p-value of 0.36. In other words, the irregularity effect was limited to the lowest values of $H_1$.

Notice however that the detection was low (61%) for fields with a Hurst anisotropy when $H_1 > 0.9$. This was not directly due to the field irregularity, but rather to the low amplitude $\delta\beta = H_2 - H_1$ of the Hurst function, which is necessarily below 0.1 when $H > 0.9$.

<table><caption>Table 5: Effect of the anisotropy amplitude on the detection.</caption><thead><tr><th colspan="2">H. anisotropy</th><th colspan="2">T. anisotropy</th></tr><tr><th>δβ</th><th>Detected cases</th><th>δτ</th><th>Detected cases</th></tr></thead><tbody><tr><td>[0, 0.01]</td><td>19/128 (15%)</td><td>[1, 1.1]</td><td>9/140 (6%)</td></tr><tr><td>[0.01, 0.02]</td><td>17/91 (19%)</td><td>[1.1, 1.2]</td><td>39/129 (30%)</td></tr><tr><td>[0.02, 0.03]</td><td>41/66 (62%)</td><td>[1.2, 1.3]</td><td>67/119 (56%)</td></tr><tr><td>[0.03, 0.04]</td><td>70/81 (86%)</td><td>[1.3, 1.4]</td><td>87/120 (72%)</td></tr><tr><td>[0.04, 0.07]</td><td>152/168 (90%)</td><td>[1.4, 1.8]</td><td>323/385 (84%)</td></tr><tr><td>[0.07, 1]</td><td>1427/1466 (97%)</td><td>[1.8, +∞]</td><td>1054/1107 (95%)</td></tr></tbody></table>
The effect of this amplitude is analyzed in more detail in Table 5. We clearly see that the detection depended on this factor up to a bound, which was about 0.07. When the amplitude was above this bound, the detection was almost perfect (97% of detected cases). The same phenomenon occurred for the detection of cases with a topothesy anisotropy. For these cases, we defined the amplitude as $δτ = C_2/C_1$. The anisotropic cases were almost perfectly detected when the amplitude was above 1.8. Hence, results obtained for these cases could have been better if we had chosen to simulate fields with higher amplitudes.

# 5 An application to microscopic images

<figure><img src="image_18.png" /><figcaption>Figure 4: Segmentation of a microscopic image. Rows 1 to 4 concern image slices at depths 5, 15, 30, and 40 µm, respectively. Columns 1 to 3 of each row show the original slice, regions with and without detected axons, respectively. Image are courtesy of Franck Debarbieux (INT, CERIMED, AMU, 2015).</figcaption></figure>
In this section, we present an application of our methodology in the context of a collaboration with the biology team of Dr. Franck Debarbieux (Institut des Neurosciences de la Timone & Centre Européen de Recherche en Imagerie Médicale, Aix-Marseille Université, France). People of this team deal with the study of inflammatory process after spinal cord injury (Fenrich et al. (2012, 2013)). For their study, they have set a protocol for imaging *in vivo* spinal cords of fluorescent transgenic mice through implanted glass windows with a two-photon microscope. Images acquired with this protocol enable one to observe some marked cells involved in inflammation at a subcellular resolution of 0.85 × 0.85 µm² and at depths varying from 5 µm to 40 µm. Figure 4 presents some image slices at different depths of a same mouse. These slices show both inflammatory cells (in a circular form) and axons (in an elongated form). The analysis of these two structures and their interactions are critical for the understanding of the inflammatory process after injury (Fenrich et al. (2012, 2013)).

So as to assist biologists in their analysis, we seek to develop automated techniques that would locate structures and characterize their interactions. The work presented here only concerns the localization part. We design a method to locate axon areas in image slices. In future work, this method is intended to serve as an initialization step for tracking axons in slices.

As illustrated in Figure 4, axons produce some anisotropic image textures whereas other structures such as cells are associated to isotropic ones. Hence, we propose to distinguish axon areas from the rest of the image by applying locally isotropy tests. More precisely, for each pixel *m* of an image, we apply an isotropy test within all square neighborhoods of size 60 × 60 containing *m*. We then compute the median of p-values of tests in all neighborhoods of *m*. If this median value is below 0.01, we consider that an axon is present in the neighborhood of *m* and label this pixel with 1. If not, the pixel is labeled with 0, meaning absence of axons. In this way, we obtain a binary segmentation of the image into axon and non-axon regions. The isotropy test is implemented with a first order-increment filter defined in Example 3.1 with $L = (0,2)$ and quadratic variations computed from scales $\sqrt{2}$ to 10 pixels (i.e. from 1.2 to 8.5 µm) in four directions (horizontal, vertical, and diagonals).

Some results are shown on Figure 4. To get images of the second (resp.
third) column, we set to zero gray-level values of the image of the first column at pixels labeled with 0 (resp. 1). We can see that major parts of axons are well separated from inflammatory cells. Misclassified axons are mainly concentrated at boundaries of axon areas. In the third and fourth slices, some cells are also misclassified due to the anisotropic aspect of their texture. We can notice that the segmentation detection is robust to intensity variations. It is also invariant to both the change of axon orientations, and the inhomogeneity of inflammatory cell textures.

# 6 Discussion

We have proposed a statistical methodology to test whether the texture of an irregular image is isotropic or not. Statistics of our isotropy tests were built upon quadratic variations. Specific to our approach, these variations were defined with increments computed in different orientations and at different scales. Hence, variations could bring useful directional information. We developed the asymptotic theory of tests in a general context of intrinsic random fields that allows us to take into account the presence of arbitrary polynomial trends in images. We established an asymptotic Gaussian linear relationship between variation and scale logarithms which involves direction-dependent intercepts. We then designed isotropy tests using Fisher statistics which check the equality of these intercepts.

The asymptotic result was obtained by applying a multivariate version of Breuer-Major (see Theorem 3.4). Proved by Biermé et al. (2011), this Breuer-Major theorem was already applied for a variation-based estimation of the Hölder irregularity in the case when the field is Gaussian with stationary increments and has a spectral density satisfying (3) with a constant Hurst function. We thus extended its application to a more general framework of IRF satisfying (3) with non-constant Hurst functions, and the analysis of multi-oriented variations.

Our test methodology is semi-parametric. It is theoretically established for a large class of Gaussian random fields whose spectral densities satisfy two weak conditions. The first is a generic non-parametric integrability condition at low frequency, the second is a condition at high frequency constraining densities to
be close to that of an AFBF. This condition involves some functional parameters (topothesy and Hurst functions). However, the test procedure does not involve any estimation of these parameters.

Some isotropy tests developed by Cabana (1987); Guan et al. (2004); Maity and Sherman (2012) for spatial data could be applied to images. These tests enable one to check if a stationary field is geometrically isotropic. However, such isotropy concerns all field frequencies whereas ours focuses on the asymptotic of the highest ones. Hence, applied to images, tests of Cabana (1987); Guan et al. (2004); Maity and Sherman (2012) would probably detect anisotropic features that are not necessarily related to textures we deal with. Moreover, they are not suitable for IRF of arbitrary orders.

We have presented an application of our statistical methodology to two-photon microscopic images. In this application, isotropic tests were used locally to segment a specific structure (axon) producing anisotropic textures. In other applications, such tests could serve as a means to detect local isotropic or anisotropic abnormalities within images (e.g. tissue pathologies, material defects, etc.). They could also be applied globally to images to classify their textures as isotropic and anisotropic.

In Image Processing, the analysis of the texture anisotropy is often associated to a detection of a privileged anisotropy direction (see Ameida (1997); Germain et al. (2003); Gorkani and Picard (1994); Lefebvre et al. (2011); Molina and Feito (2002); Rao and Schunck (1991); Kass and Witkin (1987)). This topic is beyond the scope of the work presented here, we plan to extend our methodology to include confidence intervals and tests for anisotropy directions.

# 7 Supplementary Material

In Section S1 of the supplementary material, the interested reader can find a proof of Theorem 3.4. In Section S2, propositions stated in Section 3 are also proved. Eventually, Section S3 presents the method used to estimate the covariance of quadratic variations.
Acknowledgements

We thank Dr Franck Debarbieux (Institut des Neurosciences de la Timone & Centre Européen de Recherche en Imagerie Médicale, Aix-Marseille Université, France) for the microscopic images he provided us.

References

Adler, R. J. (1981). *The Geometry of Random Field*. John Wiley & Sons.

Ameida, M. (1997). Anisotropic textures with arbitrary orientation. *J. Math. Imaging Vis.*, 7:241-251.

Aujol, J.-F., Aubert, G., Blanc-Feraud, L., and Chambolle, A. (2005). Image decomposition into a bounded variation component and an oscillating component. *J. Math. Imaging Vis.*, 22(1):71-88.

Aujol, J.-F., Gilboa, G., Chan, T., and Osher, S. (2006). Structure-texture image decomposition - modeling, algorithms and parameter selection. *Int. J. Comput. Vision*, 67(1):111-136.

Ayache, A., Bertrand, P., and Lévy-Véhel, J. (2007). A central limit theorem for the quadratic variations of the step fractional Brownian motion. *Stat. Inference Stoch. Process.*, 10:1-27.

Ayache, J. and Lévy-Véhel, J. (2004). On the identification of the pointwise Hölder exponent of the generalized multifractional Brownian motion. *Stoch. Proc. Appl.*, 111(1):119-156.

Benassi, A., Bertrand, P., Cohen, S., and Istas, J. (2000). Identification of the Hurst index of a step fractional Brownian motion. *Stat. Inference Stoch. Process.*, 3(1/2):101-111.

Benassi, A., Cohen, S., and Istas, J. (1998). Identifying the multifractional function of a Gaussian process. *Statist. Probab. Lett.*, 39:337-345.

Benson, D., Meerschaert, M. M., Bäumer, B., and Scheffler, H. P. (2006). Aquifer operator-scaling and the effect on solute mixing and dispersion. *Water Resour. Res.*, 42:1-18.

Besag, J. (1986). On the statistical analysis of dirty pictures. *J. Roy. Statist. Soc. Ser. B*, 48(3):259-302.

Biermé, H. (2005). *Champs aléatoires : autosimilarité, anisotropie et étude directionnelle*. PhD thesis, University of Orleans, France.

Biermé, H., Benhamou, C., and Richard, F. (2009). Parametric estimation for Gaussian operator scaling random fields and anisotropy analysis of bone radiograph textures. In Pohl, K., editor, *Proc. of MICCAI*, pages 13-24, London, UK.
Biermé, H., Bonami, A., and León, J. R. (2011). Central limit theorems and quadratic variations in terms of spectral density. *Electron. J. Probab.*, 16(13):362-395.

Biermé, H., Moisan, M., and Richard, F. (2015). A turning-band method for the simulation of anisotropic fractional Brownian field. *J. Comput. Graph. Statist.*, 24(3):885-904.

Biermé, H. and Richard, F. (2008). Estimation of anisotropic Gaussian fields through Radon transform. *ESAIM: Probab. Stat.*, 12(1):30-50.

Biermé, H. and Richard, F. (2011). Analysis of texture anisotropy based on some gaussian fields with spectral density. In Bergounioux, M., editor, *Mathematical Image Processing*, pages 59-73. Springer Proc.

Bonami, A. and Estrade, A. (2003). Anisotropic analysis of some Gaussian models. *J. Fourier Anal. Appl.*, 9:215-236.

Brunet-Imbault, B., Lemineur, G., Chappard, C., et al. (2005). A new anisotropy index on trabecular bone radiographic images using the fast Fourier transform. *BMC Med. Imaging*, 5:4-15.

Cabana, E. (1987). Affine processes: A test of isotropy based on level set. *SIAM J. Appl. Math.*, 47(4):886-891.

Chan, G. and Wood, T. (2000). Increment-based estimators of fractal dimension for two-dimensional surface data. *Stat. Sinica*, 10:343-376.

Chilès, J. and Delfiner, P. (2012). *Geostatistics: Modeling Spatial Uncertainty*. J. Wiley.

Coeurjolly, J. (2005). Identification of multifractional Brownian motion. *Bernoulli*, 11(6):987-1008.

Cressie, N. (1993). *Statistics for Spatial Data*. J. Wiley, New York.

Cross, G. and Jain, A. (1983). Markov random field texture models. *IEEE Trans. Pattern Anal. Mach. Intell.*, 5(1):25-39.

Davies, E. (2008). *Handbook of Texture Analysis*, chapter Introduction to Texture Analysis. World Scientific.

Davies, S. and Hall, P. (1999). Fractal analysis of surface roughness by using spatial data. *J. Roy. Statist. Soc. Ser. B*, 61:3-37.

Demanet, L. and Ying, L. (2007). Wave atoms and sparsity of oscillatory patterns. *Appl. Comput. Harmon. Anal.*, 23:368-387.

Desolneux, A., Moisan, L., and Ronsin, R. (2012). A compact representation of random phase
and Gaussian textures. In Proc. ICASSP, pages 1381–1384.

Efros, A. and Leung, T. (1999). Texture synthesis by non-parametric sampling. In Proc. ICCV, pages 1033–1038, Kerkyra, Greece.

Fenrich, K., Weber, P., Rougon, G., and Debarbieux, F. (2012). Long-term *in vivo* imaging of normal and pathological mouse spinal cord with subcellular resolution using implanted glass windows. J. Physiol., 590(16):3665–3675.

Fenrich, K., Weber, P., Rougon, G., and Debarbieux, F. (2013). Long- and short-term intravital imaging reveals differential spatiotemporal recruitment and function of myelomonocytic cells after spinal cord injury. J. Physiol., 591(19):4985–4902.

Galerne, B., Gousseau, Y., and Morel, J. (2011). Random phase textures: Theory and synthesis. IEEE Trans. Image Process., 20:257–267.

Gelfand, I. and Vilenkin, N. Y. (1964). *Generalized Functions*, volume 4: Applications to Harmonic Analysis. Academic Press.

Geman, S. and Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the bayesian restoration of images. *IEEE Trans. Pattern Anal. Mach. Intell.*, 6:721–741.

Germain, C., Da Costa, J. P., Lavialle, O., and Baylou, P. (2003). Multiscale estimation of vector field anisotropy application to texture characterization. *Signal Process.*, 83(7):1487–1503.

Gorkani, M. and Picard, R. (1994). Texture orientation for sorting photos "at a glance". In *TR-292, M.I.T., Media Laboratory, Perceptual Computing Section*, pages 459–464.

Guan, Y., Sherman, M., and Calvin, J. (2004). A nonparametric test for spatial isotropy using subsampling. *J. Am. Stat. Assoc.*, 99(467):810–821.

Haralick, R. (1979). Statistical and structural approaches to texture. Proc. IEEE, 67(5):786–804.

Istas, J. and Lang, G. (1997). Quadratic variations and estimation of the local Holder index of a Gaussian process. Ann. Inst. Henri Poincaré, Prob. Stat., 33(4):407–436.

Jain, A. and Farrokhsia, F. (1991). Unsupervised texture segmentation using gabor filters. Pattern Recogn., 24:1167–1186.

Jain, A., Hong, L., and Bolle, R. (1997). On-line fingerprint verification. IEEE Trans. Pattern Anal. Mach. Intell., 19:302–314.

Jiang, X. (2005). On orientation and anisotropy estimation for online fingerprint authentication. IEEE Trans. Signal Process., 53(10):4038–4049.
Julesz, B. and Bergen, J. (1983). Textons, the fundamental elements in preattentive vision and perception of textures. *Bell System Techn. J.*, 62(6):1619–1645.

Kaplan, L. and Kuo, J. C.-C. (1995). Texture roughness analysis and synthesis via extended self-similar (ESS) model. *IEEE Trans. Pattern Anal. Mach. Intell.*, 17(11):1043–1056.

Kass, M. and Witkin, A. (1987). Analyzing oriented patterns. *Comput. Vis. Graph. Image Process.*, 37:362–385.

Kent, J. T. and Wood, A. T. A. (1997). Estimating the fractal dimension of a locally self-similar Gaussian process by using increments. *J. Roy. Statist. Soc. Ser. B*, 59(3):679–699.

Lefebvre, A., Corpetti, T., and Moy, L. H. (2011). Estimation of the orientation of textured patterns via wavelet analysis. *Pattern Recogn. Lett.*, 32:190–196.

Lemineur, G., Harba, R., Jennane, R., et al. (2004). Fractal anisotropy measurement of bone texture radiographs. In *Proc. ISCCSP*, pages 275–278.

Maity, A. and Sherman, M. (2012). Testing for spatial isotropy under general design. *J. Stat. Plan. Infer.*, 142:1081–1091.

Mandelbrot, B. B. and Van Ness, J. (1968). Fractional Brownian motion, fractional noises and applications. *SIAM Rev.*, 10:422–437.

Matheron, G. (1973). The intrinsic random functions and their applications. *Ad. Appl. Prob.*, 5:439–468.

Meyer, Y. (2001). Oscillating patterns in image processing and in some nonlinear evolution equation. In *The 15th Dean Jacqueline B. Lewis Memorial Lectures*.

Molina, A. and Feito, F. (2002). A method for testing anisotropy and quantifying its direction in digital images. *Comput. Graphics*, 26(771–784).

Pentland, A. (1984). Fractal-based description of natural scenes. *IEEE Trans. Pattern Anal. Mach. Intell.*, 6:661–674.

Peyré, G. (2009). Manifold models for signals and images. *Comput. Vis. Image Und.*, 113(2):249–260.

Portilla, J. and Simoncelli, E. (2000). A parametric texture model based on joint statistics of complex wavelet coefficients. *Int. J. Comput. Vision*, 40(1):49–71.

Rao, R. and Lohse, G. (1993). Identifying high-level features of texture perception. *Comput. Vis. Graph. Image Process.*, 55(2):218–233.

Rao, R. and Schunck, B. (1991). Computing oriented texture fields. In *Comput. Vis. Graph*.
*Image Process.*, volume 53, pages 157–185.

Richard, F. (2015). Analysis of anisotropic Brownian textures and application to lesion detection in mammograms. *Procedia Environ. Sci.*, 27(16-20).

Richard, F. and Biermé, H. (2010). Statistical tests of anisotropy for fractional Brownian textures. Application to full-field digital mammography. *J. Math. Imaging Vis.*, 36(3):227–240.

Tuceryan, M. and Jain, A. (1998). *The handbook of Pattern Recognition and Computer Vision*, chapter Texture analysis. World Scientific.

Unser, M. (1995). Texture segmentation and classification using wavelet frames. *IEEE Trans. Image Process.*, 4(11):1549–1560.

Wei, L.-Y., Lefebvre, S., Kwatra, V., and Turk, G. (2009). State of the art in example-based texture synthesis. In *Eurographics 2009*, State of the Art Report. Eurographics Association, EG-STAR.

Yaglom, A. (1986). *Correlation Theory of Stationary and Related Random Functions: Basics Results*. Springer Verlag.

Zhu, S., Wu, Y., and Mumford, D. (1998). Filters, random fields and maximum entropy (frame): Towards a unified theory for texture modeling. *Int. J. Comput. Vision*, 27(2):107–126.

Zhu, Z. and Stein, M. (2002). Parameter estimation for fractional Brownian surfaces. *Stat. Sinica*, 12:863–883.

Author affiliation: Aix Marseille Université, CNRS, Centrale Marseille, I2M, UMR 7373, 13453
Marseille, France. E-mail: frederic.richard@univ-amu.fr