Patch redundancy in images: a statistical testing framework and some applications

TL;DR

This paper introduces a statistical testing framework to analyze local spatial redundancy in natural images, enabling applications like denoising, periodicity detection, and texture ranking through a fast algorithm.

Contribution

The work develops a novel a contrario statistical model for patch similarity, providing a rigorous criterion for redundancy detection in images.

Findings

Effective redundancy detection algorithm

Applications in denoising and texture analysis

Non-asymptotic probability expressions for similarity measures

Abstract

In this work we introduce a statistical framework in order to analyze the spatial redundancy in natural images. This notion of spatial redundancy must be defined locally and thus we give some examples of functions (auto-similarity and template similarity) which, given one or two images, computes a similarity measurement between patches. Two patches are said to be similar if the similarity measurement is small enough. To derive a criterion for taking a decision on the similarity between two patches we present an a contrario model. Namely, two patches are said to be similar if the associated similarity measurement is unlikely to happen in a background model. Choosing Gaussian random fields as background models we derive non-asymptotic expressions for the probability distribution function of similarity measurements. We introduce a fast algorithm in order to assess redundancy in natural images and present applications in denoising, periodicity analysis and texture ranking.

Figures

Equations

$$\mathcal{A}\mathcal{S}(u,\mathrm{t},\omega )=\Vert P_{\mathrm{t}+\omega }(u)-P_{\omega }(u){\Vert }_2^2.$$

$$\mathsf{AP}(\mathrm{t},\omega ,a)={\mathbb{P}}_0\left[\mathcal{A}\mathcal{S}(U,\mathrm{t},\omega )\le a(\mathrm{t})\right].$$

$$\mathsf{ANFA}(\omega ,a)=\sum _{\mathrm{t}\in \Omega }\mathsf{AP}(\mathrm{t},\omega ,a).$$

$$\mathsf{AP}\left(\mathrm{t},\omega ,{\mathsf{AP}}^{-1}\left(\mathrm{t},\omega ,q\right)\right)=q.$$

$$\mathsf{ANFA}(\omega ,a)={\mathrm{NFA}}_{\text{max}}\text{and}{\mathbb{P}}_0\left[\text{“at least n offsets are detected in U” }\right]\le \frac{{\mathrm{NFA}}_{\text{max}}}{n}.$$

$$\mathsf{ANFA}(\omega ,a)=\sum _{\mathrm{t}\in \Omega }\mathsf{AP}(\mathrm{t},\omega ,a)=\sum _{\mathrm{t}\in \Omega }\mathsf{AP}\left(\mathrm{t},\omega ,{\mathsf{AP}}^{-1}\left(\mathrm{t},\omega ,{\mathrm{NFA}}_{\text{max}}/|\Omega |\right)\right)={\mathrm{NFA}}_{\text{max}},$$

$${\mathbb{P}}_0\left[\text{“at least n offsets are detected in U” }\right]$$

$$\le \frac{{\sum }_{\mathrm{t}\in \Omega }\mathbb{E}\left[1_{\mathcal{A}\mathcal{S}(U,\mathrm{t},\omega )\le a(\mathrm{t})}\right]}{n}\le \frac{{\mathrm{NFA}}_{\text{max}}}{n},$$

$$\mathcal{A}\mathcal{S}(u,\mathrm{t},\omega )\le {\mathsf{AP}}^{-1}\left(\mathrm{t},\omega ,{\mathrm{NFA}}_{\text{max}}/|\Omega |\right).$$

$${\mathbb{P}}_0\left[\mathcal{A}\mathcal{S}(U,\mathrm{t},\omega )\le \mathcal{A}\mathcal{S}(u,\mathrm{t},\omega )\right]=\mathsf{AP}\left(\mathrm{t},\omega ,\mathcal{A}\mathcal{S}(u,\mathrm{t},\omega )\right)\le {\mathrm{NFA}}_{\text{max}}/|\Omega |.$$

$$m_u=\sum _{\mathrm{x}\in \Omega }u(x)/|\Omega |,\text{and}U=|\Omega {|}^{-1/2}(u-m_u)\ast W.$$

$$\mathbb{E}\left[U(\mathrm{x})\right]=0\text{and}\mathrm{Cov}\left[U(\mathrm{x}),U(\mathrm{y})\right]=|\Omega {|}^{-1}\sum _{\mathrm{z}\in \Omega }(\dot{u}(\mathrm{z})-m_u)(\dot{u}(\mathrm{z}-(\mathrm{y}-\mathrm{x}))-m_u).$$

$${\Delta }_f(\mathrm{t},\mathrm{x})=2{\Gamma }_f(\mathrm{x})-{\Gamma }_f(\mathrm{x}+\mathrm{t})-{\Gamma }_f(\mathrm{x}-\mathrm{t}).$$

$$T=\left\{\mathrm{t}\in {\mathbb{Z}}^2,\mathrm{t}+\omega \subset \Omega ,\Vert \mathrm{t}{\Vert }_{\infty }\le c\right\},$$

$$\hat{u}(\mathrm{x})={\left|\{\mathrm{t}\in \Omega ,\text{s.t}\mathrm{x}\in \mathrm{t}+\omega \subset \Omega \}\right|}^{-1}\sum _{\mathrm{t}\in \Omega ,\text{s.t}\mathrm{x}\in \mathrm{t}+\omega \subset \Omega }\hat{p}(u,\mathrm{t}+\omega )(\mathrm{x}).$$

$$\hat{p}(u,\omega )=\sum _{\mathrm{t}\in T}{\lambda }_{\mathrm{t}}{P}_{\mathrm{t}+\omega }(u),{\lambda }_{\mathrm{t}}=\frac{1_{\mathcal{A}\mathcal{S}(u,\mathrm{t},\omega )\le a(\mathrm{t})}}{{\sum }_{\mathrm{s}\in T}{1}_{\mathcal{A}\mathcal{S}(u,\mathrm{s},\omega )\le a(\mathrm{s})}}.$$

$${\lambda }_{\mathrm{t}}=\frac{\exp \left(-\frac{\mathcal{A}\mathcal{S}(u,\mathrm{t},\omega )}{h^2}\right)}{{\sum }_{\mathrm{t}\in T}\exp \left(-\frac{\mathcal{A}\mathcal{S}(u,\mathrm{t},\omega )}{h^2}\right)}.$$

$$a(\mathrm{t})={\mathsf{AP}}^{-1}\left(\mathrm{t},\omega ,1-{\mathrm{NFA}}_{\text{max}}/|T|\right),$$

$${\mathbb{P}}_0\left[|T|-N_{\omega }(W)\ge n\right]\le \frac{{\mathrm{NFA}}_{\text{max}}}{n}.$$

$${\mathbb{P}}_0\left[|T|-N_{\omega }(W)\ge n\right]\le \frac{|T|-{\sum }_{\mathrm{t}\in T}\mathbb{E}\left[1_{\mathcal{A}\mathcal{S}(W,\mathrm{t},\omega )\le a(\mathrm{t})}\right]}{n}\le \frac{{\mathrm{NFA}}_{\text{max}}}{n}.$$

$$\mathbb{P}\left[\Vert \hat{p}(U,\omega )-P_{\omega }(u_0){\Vert }_2\le \sigma (a_T^{1/2}+a_W^{1/2})|\hat{T}\right]\ge 1-{\epsilon }_W.$$

$$\Vert P_{\mathrm{t}+\omega }(U)-P_{\omega }(u_0){\Vert }_2$$

$$\le \Vert P_{\mathrm{t}+\omega }(U)-P_{\omega }(U){\Vert }_2+\Vert P_{\omega }(U)-P_{\omega }(u_0){\Vert }_2$$

$$\le \sigma a_T^{1/2}+\sigma \Vert P_{\omega }(W){\Vert }_2.$$

$$\left\{\Vert P_{\omega }(W){\Vert }_2\le a_W^{1/2}\right\}\subset \left\{\Vert P_{\mathrm{t}+\omega }(U)-P_{\omega }(u_0){\Vert }_2\le \sigma (a_T^{1/2}+a_W^{1/2})\right\},$$

$$\mathbb{P}\left[\Vert \hat{p}(U,\omega )-P_{\omega }(u_0){\Vert }_2\le \sigma (a_T^{1/2}+a_W^{1/2})|\hat{T}\right]$$

$$\ge \mathbb{P}\left[\bigcap _{\mathrm{t}\in \hat{T}}\{\Vert P_{\mathrm{t}+\omega }(U)-P_{\omega }(u_0){\Vert }_2^2\le {\sigma }^2(a_T^{1/2}+a_W^{1/2}{)}^2\}|\hat{T}\right]$$

$$\ge \mathbb{P}\left[\Vert P_{\omega }(W){\Vert }_2^2\le a_W|\hat{T}\right]\ge 1-{\epsilon }_W.$$

$$C_{\mathrm{t}}=\left(\begin{array}{cccc}{B}_0& B_1& \dots & B_{p-1}\\ B_1^{\top }& B_0& \ddots & ⋮\\ ⋮& \ddots & B_0& B_1\\ B_{p-1}^{\top }& \dots & B_1^{\top }& B_0\end{array}\right)+2\mathrm{Id},\{\begin{array}{ll}{B}_{\ell }=D_{|t_y|}\in {\mathcal{M}}_p(\mathbb{R})& \text{if }\ell =|t_x|\\ B_{\ell }=0& \text{otherwise}\end{array}$$

$${\mathcal{I}}_{\omega }(\mathrm{t})$$

$$=\sum _{\mathrm{x}\in \omega }(\dot{u}(\mathrm{x})-\dot{u}(\mathrm{x}+\mathrm{t}){)}^2=\mathcal{A}\mathcal{S}(u,\mathrm{t},\omega ).$$

$$\mathcal{L}(B,M,{\sigma }^2|E)=-2(|E|+1)\log ({\sigma }^2)-\frac{1}{2{\sigma }^2}\underset{q(B,M|E)}{\underbrace{\left(\sum _{\mathrm{e}\in E}\Vert m_{\mathrm{e}}{b}_1+n_{\mathrm{e}}{b}_2-\mathrm{e}{\Vert }^2+r(B,M)\right)}},$$

$$\tilde{M}={\left({\Lambda }_{B_n}\otimes {\mathrm{Id}}_{|E|}\right)}^{-1}{E}_{B_n}\in {\mathbb{R}}^{2|E|},B_{n+1}={\left({\Lambda }_{M_{n+1}}\otimes {\mathrm{Id}}_2\right)}^{-1}{E}_{M_{n+1}}\in {\mathbb{R}}^4,$$

Full text

K****eywords patch, redundancy, statistical framework, a contrario method, image denoising, texture, periodicity analysis.

1 Introduction

In many image processing applications, using local information combined with the knowledge of long-range spatial arrangement is crucial. The spatial redundancy on sub-images called patches, encodes the small scale structure of the image as well as its large scale organization. More precisely, local information is encoded in the patch content and the large scale organization is contained in the redundancy of this information across the patches of the image. For example, patch-based inpainting techniques, such as [10, 33], assign patches of a known region to patches of an unknown region. Namely, each patch position on the border of the unknown region is associated to an offset corresponding to the best patch according to the partial available information. In [33] the authors replace the search on the whole image by a search among the most redundant offsets in the known region. This allows the authors of [33] to retrieve long-range spatial structure in the unknown part of the image. Another famous application of spatial redundancy can be found in denoising, with the seminal work (Non-Local means) of Buades and coauthors [5], in which the authors propose to replace a noisy patch by the mean over all spatially redundant patches.

Last but not least, spatial redundancy is of crucial importance in exemplar-based texture synthesis. In this paper we define textures as images containing repeated patterns but also reflecting randomness in the arrangement of these patterns. Among textures, one important class is given by the microtextures in which no individual object can be clearly delimited. In the periodic case, a more precise definition will be given in Definition 4. These microtexture models can be described by Gaussian random fields [62, 27, 42, 68]. Parametric models using features such as wavelet transform coefficients [55], scattering transform coefficients [59] or convolutional neural network outputs [29] have been proposed in order to derive image models with more structure. On the other hand, non-parametric patch-based algorithms such as [25, 24, 38, 56, 28] propose to use most similar patches in order to fill the new texture images, similarly to inpainting techniques.

All these techniques lift images in spaces with dimensions higher than the original image space, and make use of the redundancy of the lifting to extract important structural information. There exist two main types of lifting: feature extraction or patch extraction. Feature extraction relies on the use of filters, linear or non-linear, which aim at selecting substantial local information. Among popular kernels are oriented and multiscale filters, which happened to be identified as early processing in mammal vision systems [13, 35]. These last years have seen the rise of neural networks in which the filter dictionary is no longer given as an input but learned through a data-driven optimization procedure [60]. On the other hand, patch-based methods rely on the assumption that image processing tasks are simplified when conducted in the higher dimensional patch space.

Every analysis performed in a lifted space, built via feature extraction or patch extraction, relies on the comparison of points in this space. In patch-based lifted spaces, we aim at finding dissimilarity functions such that two patches are visually close if the dissimilarity measurement between them is small. In this paper we focus on the square Euclidean distance but other choices could be considered [64, 65, 15, 17].

This leads us to consider a statistical hypothesis testing framework to assess similarity (or dissimilarity) between patches. The null hypothesis is defined as the absence of local structural similarities in the image. Reciprocally the alternative hypothesis is defined as the presence of such similarities. There exists a wide variety of tractable models exhibiting no similarity at long-range, like Gaussian random fields [62, 27, 42, 68] or spatial Markov random fields [11], whereas sampling and inference in very structured models rely on optimization procedures and may be computationally expensive, their distribution being the limit of some Markov chain [70, 47] or some stochastic optimization procedure [4]. This encourages us to consider an a contrario approach, i.e. we do not consider the alternative hypothesis and focus on rejecting the null hypothesis. This framework was successfully applied in many areas of image processing [14, 19, 20, 1, 7] and aims at identifying structure events in images. This statistical model takes its roots in the fundamental work of the Gestalt theory [21]. One of its principle, the non-accidentalness principle [46] or Helmholtz principle [69, 20], states that no structure is perceived in a noise model. To be precise, in our case of interest, we want to assess that no spatial redundancy is perceived in microtexture models. This methodology allows us to only design a locally structured background model to define a null hypothesis. Combining a contrario principles and patch-based measures, we propose an algorithm to identify auto-similarities in images.

We then turn to the implementation of such an algorithm and illustrate the diversity of its possible applications with three examples: denoising, lattice extraction, and periodicity ranking of textures. In our denoising application we propose a modification of the celebrated Non-Local means algorithm [5] (NL-means) by inserting a threshold in the selection of similar patches. Using an a contrario model we are able to give probabilistic control on the patch reconstruction.

We then focus on periodicity detection and, more precisely, lattice extraction. Periodicity in images was described as an important feature in early mathematical vision [32]. Most of the proposed methods to analyze periodicity rely on global measurements such as the modulus of the Fourier transform [49] or the autocorrelation [43]. These global techniques are widely used in crystallography where lattice properties, such as the angle between basis vectors, are fundamental [50, 58]. Since all of our measurements are local, we are able to identify periodic similarities even in images which are not periodic but present periodic parts, for instance if two crystal structures are present in a single crystallography image. We draw a link between the introduced notion of auto-similarity and the inertia measurement in co-occurence matrices [32]. We then introduce our lattice proposal algorithm which combines a detection map, i.e. the output of our redundancy detection algorithm, and graphical model techniques, as in [53], in order to extract lattice basis vectors.

Our last application concerns texture ranking. Since the definition of texture is broad and covers a wide range of images, it is a natural question to identify criteria in order to distinguish textures. In [45], the authors use a classical measure for distinguishing textures: regularity. In this work, we narrow this criterion and restrict ourselves to the study of periodicity in texture images. The proposed graphical model inference naturally gives a quantitative measurement for texture periodicity ranking. We give an example of ranking on 25 images of the Brodatz set.

Our paper is organized as follows. An a contrario framework for local similarity detection is proposed in Section 1. In the a contrario framework, a background model, corresponding to the null hypothesis, is required. The consequence of choosing Gaussian models as background models is investigated and a redundancy detection algorithm is proposed in Section 3. The rest of the paper is dedicated to some examples of application of the proposed framework. After reviewing one of the most popular method in image denoising we introduce a denoising algorithm in Section 4.1 and present our experimental results in Section 4.2. Local dissimilarity measurements can be used as periodicity detectors. The link between the locality of the introduced functions and the literature on periodicity detection problems is investigated in Section 5.1. An algorithm for detecting lattices in images is given in Section 5.2 and numerical results are presented in Section 5.3. In our last experiment in Section 5.4, we introduce a criterion for measuring texture periodicity. We conclude our study and discuss future work in Section 6.

2 An a contrario framework for auto-similarity

We first introduce a notion of dissimilarity between patches of an input image.

Definition 1** (Auto-similarity)**

Let $u$ be an image defined over a domain $\Omega=\llbracket 0,M-1\rrbracket^{2}\subset\mathbb{Z}^{2}$, with $M\in\mathbb{N}\backslash\{0\}$. Let $\omega\subset\mathbb{Z}^{2}$ be a patch domain. We introduce ${P_{\omega}(u)=(\dot{u}(\bm{\mathrm{y}}))_{\bm{\mathrm{y}}\in\omega}}$ the patch at position $\omega$ in the periodic extension of $u$ to $\mathbb{Z}^{2}$, denoted by $\dot{u}$. We define the auto-similarity with patch domain $\omega$ and offset $\bm{\mathrm{t}}\in\mathbb{Z}^{2}$ by

[TABLE]

The auto-similarity computes the distance between a patch of $u$ defined on a domain $\omega$ and the patch of $u$ defined by the domain $\omega$ shifted by the offset vector $\bm{\mathrm{t}}$.

In what follows, we introduce an a contrario framework on the auto-similarity. This framework will allow us to derive an algorithm for detecting spatial redundancy in natural images.

In this section we fix an image domain $\Omega\subset\mathbb{Z}^{2}$ and a patch domain $\omega\subset\Omega$. We recall that our final aim is to design a criterion that will answer the following question: are two given patches similar? This criterion will be given by the comparison between the value of a dissimilarity function and a threshold $a$. We will define the threshold $a$ so that few similarities are identified in the null hypothesis model, i.e. similarity does not occur “just by chance”. Thus we can reformulate the initial question: is the similarity output of a dissimilarity function between two patches small enough? Or, to be more precise, how can we set the threshold $a$ in order to obtain a criterion for assessing similarity between patches?

This formulation agrees with the a contrario framework [21] which states that geometrical and/or perceptual structure in an image is meaningful if it is a rare event in a background model. This general principle is sometimes called the Helmholtz principle [69] or the non-accidentalness principle [46]. Therefore, in order to control the number of similarities identified in the background model, we study the probability density function of the auto-similarity function with input random image $U$ over $\Omega$. We will denote by $\mathbb{P}_{0}$ the probability distribution of $U$ over $\mathbb{R}^{\Omega}$, the images over $\Omega$. We will assume that $\mathbb{P}_{0}$ is a microtexture model, see Definition 4 below for a precise definition of such a model. We define the following significant event which encodes spatial redundancy: $\mathcal{AS}(u,\bm{\mathrm{t}},\omega)\leq a(\bm{\mathrm{t}})$, where $a$, the threshold function, is defined over the offsets ($\bm{\mathrm{t}}\in\mathbb{Z}^{2}$) but also depends on other parameters such as $\omega$ or $\mathbb{P}_{0}$ . The dependency of $a$ with respect to $\bm{\mathrm{t}}$ cannot be omitted. For instance, even in a Gaussian white noise $W$, the probability distribution function of $\mathcal{AS}(W,\bm{\mathrm{t}},\omega)$ depends on $\bm{\mathrm{t}}$.

The Number of False Alarms ($\operatorname{NFA}$) is a crucial quantity in the a contrario methodology. A false alarm is defined as an occurrence of the significant event in the background model $\mathbb{P}_{0}$. We recall that in our model the significant event is patch redundancy. This test must be conducted for every possible configurations of the significant event, i.e. in our case we test every possible offset $\bm{\mathrm{t}}$. The $\operatorname{NFA}$ is then defined as the expectation of the number of false alarms over all possible configurations. Bounding the $\operatorname{NFA}$ ensures that the probability of identifying $k$ offsets with spatial redundancy is also bounded, see Proposition 1. In what follows we give the definition of the $\operatorname{NFA}$ in the spatial redundancy context.

Definition 2** ($\operatorname{NFA}$)**

Let $U\sim\mathbb{P}_{0}$, where $\mathbb{P}_{0}$ is a background microtexture model. We define the auto-similarity probability map $\mathsf{AP}$ for any $\bm{\mathrm{t}}\in\Omega$, $\omega\subset\Omega$ and $a\in\mathbb{R}^{\Omega}$ by

[TABLE]

We define the auto-similarity expected number of false alarms $\mathsf{ANFA}$ by

[TABLE]

Note that $\mathsf{AP}(\bm{\mathrm{t}},\omega,a)$ corresponds to the probability that $\omega+\bm{\mathrm{t}}$ is similar to $\omega$ in the background model $U$. For any $\bm{\mathrm{t}}\in\Omega$, the cumulative distribution function of the auto-similarity random variable $\mathcal{AS}(U,\bm{\mathrm{t}},\omega)$ under $\mathbb{P}_{0}$ evaluated at value $\alpha(\bm{\mathrm{t}})$ is given by $\mathsf{AP}(\bm{\mathrm{t}},\omega,\alpha(\bm{\mathrm{t}}))$. We denote by ${q\mapsto\mathsf{AP}^{-1}(\bm{\mathrm{t}},\omega,q)}$ the inverse cumulative distribution function, potentially defined by a generalized inverse ($\mathsf{AP}^{-1}(\bm{\mathrm{t}},\omega,q)=\inf\{\alpha(\bm{\mathrm{t}})\in\mathbb{R},\ \mathsf{AP}(\bm{\mathrm{t}},\omega,\alpha(\bm{\mathrm{t}}))\geq q\}$), of the auto-similarity random variable for a fixed offset $\bm{\mathrm{t}}$, with $q\in(0,1)$ a quantile. We now have all the tools to control the number of detected offsets in the background model.

Definition 3** (Detected offset)**

Let $u\in\mathbb{R}^{\Omega}$ be an image, $\omega\subset\Omega$ a patch domain, and $a\in\mathbb{R}^{\Omega}$. An offset $\bm{\mathrm{t}}$ is said to be detected with respect to $a$, if $\mathcal{AS}(u,\bm{\mathrm{t}},\omega)\leq a(\bm{\mathrm{t}})$.

Note that a detected offset in $U\sim\mathbb{P}_{0}$ corresponds to a false alarm in the a contrario model. In what follows we suppose that the cumulative distribution function of $\mathcal{AS}(U,\bm{\mathrm{t}},\omega)$ is invertible for every $\bm{\mathrm{t}}\in\Omega$. This ensures that for any $\bm{\mathrm{t}}\in\Omega$ and $q\in(0,1)$ we have

[TABLE]

Proposition 1

Let $\operatorname{NFA}_{\text{max}}\geq 0$ and for all $\bm{\mathrm{t}}\in\Omega$ define $a(\bm{\mathrm{t}})=\mathsf{AP}^{-1}\left(\bm{\mathrm{t}},\omega,\operatorname{NFA}_{\text{max}}/|\Omega|\right)$. We have that for any $n\in\mathbb{N}\backslash\{0\}$,

[TABLE]

Proof:

Using (3), and $a(\bm{\mathrm{t}})=\mathsf{AP}^{-1}\left(\bm{\mathrm{t}},\omega,\operatorname{NFA}_{\text{max}}/|\Omega|\right)$, we get

[TABLE]

where the last equality is obtained using (4). Concerning the upper-bound, we have, using the Markov inequality and (2), for any $n\in\mathbb{N}\backslash\{0\}$

[TABLE]

where $\mathbb{1}_{\mathcal{AS}(U,\bm{\mathrm{t}},\omega)\leq a(\bm{\mathrm{t}})}=1$ if $\mathcal{AS}(U,\bm{\mathrm{t}},\omega)\leq a(\bm{\mathrm{t}})$ and [math] otherwise. $\square$

Thus, setting $a$ as in Proposition 1, we have that an offset $\bm{\mathrm{t}}\in\Omega$ is detected for an image $u\in\mathbb{R}^{\Omega}$ if

[TABLE]

This a contrario detection framework can then be simply rewritten as 1) computing the auto-similarity function with input image $u$, 2) thresholding the obtained dissimilarity map with the inverse cumulative distribution function of the computed dissimilarity function under $\mathbb{P}_{0}$. The computed threshold depends on the offset and Proposition 1 ensures probabilistic guarantees on the expected number of detections under $\mathbb{P}_{0}$. Using the inverse property of the inverse cumulative distribution function and (5), we obtain that an offset is detected if and only if

[TABLE]

Therefore, the thresholding operation can be conducted either on $\mathcal{AS}(u,\bm{\mathrm{t}},\omega)$, see (5), or on $\mathsf{AP}\left(\bm{\mathrm{t}},\omega,\mathcal{AS}(u,\bm{\mathrm{t}},\omega)\right)$, see (6). This property will be used in Section 3.2 to define a similarity detection algorithm based on the evaluation of $\mathcal{AS}(u,\bm{\mathrm{t}},\omega)$.

3 Gaussian model and detection algorithm

3.1 Choice of background model

In this section we compute $\mathsf{AP}\left(\bm{\mathrm{t}},\omega,\alpha\right)$, i.e. the cumulative distribution function of the similarity function under the null hypothesis model, with a Gaussian background model. Indeed, if the background model is simply a Gaussian white noise the similarities identified by the a contrario algorithm are the ones that are not likely to be present in the Gaussian white noise image model. More generally, we consider stationary Gaussian random fields defined in the following way: we introduce an image $f$ over $\mathbb{R}^{\Omega}$ which contains the microtexture information we want to discard in our a contrario model. In what follows we give the definition of the microtexture model associated to $f$.

Definition 4** (Microtexture model)**

Let $f\in\mathbb{R}^{\Omega}$, we define the associated microtexture model $U$ by setting, $U=f*W$, where $*$ is the periodic convolution operator over $\Omega$ given by $v*w(\bm{\mathrm{x}})=\sum_{\bm{\mathrm{y}}\in\Omega}\dot{v}(\bm{\mathrm{y}})\dot{w}(\bm{\mathrm{x}}-\bm{\mathrm{y}})$ and $W$ is a white noise over $\Omega$, i.e. $(W(\bm{\mathrm{x}}))_{\bm{\mathrm{x}}\in\Omega}$ are i.i.d. $\mathcal{N}(0,1)$ random variables.

Given an image $u\in\mathbb{R}^{\Omega}$, a microtexture model can be derived considering

[TABLE]

Note that if $U$ is given by (7) we have for any $\bm{\mathrm{x}},\bm{\mathrm{y}}\in\Omega$

[TABLE]

We refer to [27] for a mathematical study of this model.

3.2 Detection algorithm

In this section, $\Omega$ is a finite square domain in $\mathbb{Z}^{2}$. We fix $\omega\subset\Omega$. We also define $f$, a function over $\Omega$. We consider the Gaussian random field $U=f*W$, where $W$ is a Gaussian white noise over $\Omega$. We denote by $\Gamma_{f}$ the autocorrelation of $f$, i.e. $\Gamma_{f}=f*\check{f}$ where for any $\bm{\mathrm{x}}\in\Omega$, $\check{f}(\bm{\mathrm{x}})=f(-\bm{\mathrm{x}})$. We introduce the offset correlation function $\Delta_{f}$ defined for any $\bm{\mathrm{t}},\bm{\mathrm{x}}\in\Omega$ by

[TABLE]

The following proposition, proved in [15], gives the explicit probability distribution function of the squared $\ell^{2}$ auto-similarity.

Proposition 2** (Squared $\ell^{2}$ auto-similarity function exact probability distribution function)**

Let $\Omega=\llbracket 0,M-1\rrbracket^{2}$ with $M\in\mathbb{N}\backslash\{0\}$, $\omega\subset\Omega$, $f\in\mathbb{R}^{\Omega}$ and $U=f*W$ where $W$ is a Gaussian white noise over $\Omega$. Then, for any $\bm{\mathrm{t}}\in\Omega$, $\mathcal{AS}(U,\bm{\mathrm{t}},\omega)$ has the same distribution as $\sum_{k=0}^{|\omega|-1}{\lambda_{k}(\bm{\mathrm{t}},\omega)Z_{k}}$, with $Z_{k}$ independent chi-square random variables with parameter 1 and $\lambda_{k}(\bm{\mathrm{t}},\omega)$ the eigenvalues of the covariance matrix $C_{\bm{\mathrm{t}}}$ associated with function $\Delta_{f}(\bm{\mathrm{t}},\cdot)$ restricted to $\omega$, defined in (9), i.e  for any $\bm{\mathrm{x_{1}}},\bm{\mathrm{x_{2}}}\in\omega$, $C_{\bm{\mathrm{t}}}(\bm{\mathrm{x_{1},x_{2}}})=\Delta_{f}(\bm{\mathrm{t}},\bm{\mathrm{x_{1}-x_{2}}})$.

As a consequence if $f=\delta_{0}$, i.e. $U$ is a Gaussian white noise, and $\{\bm{\mathrm{x}}+\bm{\mathrm{t}},\bm{\mathrm{x}}\in\omega\}\cap\omega=\emptyset$, i.e. there is no overlapping between the patch domain $\omega$ and its shifted version, then $\mathcal{AS}(U,\bm{\mathrm{t}},\omega)$ is a chi-square random variable with parameter $|\omega|$.

In order to compute the cumulative distribution function of a quadratic form of Gaussian random variables we must deal with two issues: 1) the computation of the eigenvalues $\lambda_{k}(\bm{\mathrm{t}},\omega)$ might be time-consuming and efficient methods must be developed ; 2) the exact computation of the cumulative distribution function of a quadratic form of Gaussian random variables requires the use of heavy integrals, see [36]. In [15] a projection method is introduced in order to easily compute approximated eigenvalues, with equality when $\omega=\Omega$. The so-called Wood F method (see [66, 3]) shows the best trade-off between accuracy and computational cost to approximate the cumulative distribution function of quadratic forms in Gaussian random variables with given weights. It is a moment method of order 3, fitting a Fisher-Snedecor distribution to the empirical one. Note that in [44] another moment method of order 3 is proposed. In what follows, we assume that we can compute the cumulative distribution function of $\mathcal{AS}(U,\bm{\mathrm{t}},\omega)$ and we refer to [15] for further details.

In Algorithm 1 we propose an a contrario framework for spatial redundancy detection. We suppose that $u$ and $\omega$ are provided by the user. Using Proposition 1 and (6) , we say that an offset is detected if $\mathsf{AP}\left(\bm{\mathrm{t}},\omega,\mathcal{AS}(u,\bm{\mathrm{t}},\omega)\right)\leq\operatorname{NFA}_{\text{max}}/|\Omega|$. The value $\operatorname{NFA}_{\text{max}}$ is supposed to be set by the user. The background model used in the auto-similarity detection is the one given in (7). Therefore, Proposition 2 and the discussion that follows can be used to compute an approximation of $\mathsf{AP}(\bm{\mathrm{t}},\omega,\mathcal{AS}(u,\bm{\mathrm{t}},\omega))$. In Figure 2 we apply Algorithm 1 to a texture image.

4 Denoising

4.1 NL-means and a contrario framework

In this section we apply the a contrario framework to the context of image denoising and propose a simple modification of the celebrated image denoising algorithm Non-Local Means (NL-means). This algorithm was introduced in the seminal paper of Buades et al. [5] and was inspired by the work of Efros and Leung in texture synthesis [25]. It was also independently introduced in [2]. This algorithm relies on the simple idea that denoising operations can be conducted in the lifted patch space. In this space the usual Euclidean distance acts as a good similarity detector and we can obtain a denoised patch by averaging all the patches with weights that depend on this Euclidean distance. Usually the weight function is set to have exponential decay, but it was suggested in [30, 57, 22] to use compactly supported weight functions in order to avoid the loss of isolated details. Since its introduction, many algorithms derived from NL-means have been proposed in order to embed the algorithm in general statistical frameworks [23, 41] or to take into account the underlying geometry of the patch space [34]. Among the state-of-the-art denoising algorithms, see [39] for a review, we consider Block-Matching and 3D Filtering (BM3D) [12] to compare our algorithm with.

There exist several works combining a contrario models and denoising tasks. Coupier et al. in [9] propose to combine morphological filters and a testing hypothesis framework to remove impulse noise. In [18] Delon and Desolneux compare different statistical frameworks to perform denoising with Gaussian noise or impulse noise. The a contrario model was also successfully used to deal with speckle noise [26] and quasi-periodic noise [61], and rely on the thresholding of wavelet or Fourier coefficients. In [37], Kervrann and Boulanger derive approximated probabilistic thresholds using $\chi_{2}$ probability distribution functions. In [67] the authors propose a testing framework in order to estimate thresholds. The expressions they derive also relies on an approximation of the probability distribution of the squared Euclidean norm between two patches in Gaussian white noise.

Following a standard extension procedure of the NL-means algorithm we consider a threshold version of it, see Algorithm 2. In what follows we fix a “clean”, or original, image $u_{0}$ defined over $\Omega$, a finite rectangular domain of $\mathbb{Z}^{2}$, a noisy image $u=u_{0}+\sigma w$, with $w$ a realization of a standard Gaussian random field $W$ and $\sigma>0$ the standard deviation of the noise. In all of our experiments we suppose that $\sigma$ is known. Note that there exist several algorithms to estimate $\sigma$ from real images, see [54] for instance. Our goal is to retrieve $u_{0}$ based on the information in $u$. We consider the lifted version of $u$ in a patch space. Let $\omega_{0}$ be a centered $8\times 8$ patch domain. For a patch window $\omega=\bm{\mathrm{x}}+\omega_{0}$ the patch search window $T$ will be defined by

[TABLE]

with $c\in\mathbb{N}$. $|T|$ denotes the cardinality of $T$. There exists a large literature concerning the setting of $c$ and $\omega_{0}$, see [22]. Note that the locality of the patch window was assessed to be a crucial feature of NL-means [31]. Suppose we have a collection of denoised patches $\hat{p}(u,\omega)$ for all patch domains $\omega$, we obtain a pixel at position $\bm{\mathrm{x}}$ in the denoised image $\hat{u}$ using the following average, see [6],

[TABLE]

We now introduce our modification of NL-means. We suppose that we are provided a threshold function $a$. The choice of such a function is discussed in Proposition 3.

Note here that the output denoised version of the patch $\hat{p}(u,\omega)$ verifies the following equation

[TABLE]

In the original NL-means method, we have

[TABLE]

Setting $h$ is not trivial and depends on many parameters (patch size, search window size, content of the original image). As in Algorithm 2, we denote $N_{\omega}(u)=\sum_{\bm{\mathrm{t}}\in T}\mathbb{1}_{\mathcal{AS}(u,\bm{\mathrm{t}},\omega)\leq a(\bm{\mathrm{t}})}$. The following proposition, similar to Proposition 1, gives a method for setting $a$. We say that an offset $\bm{\mathrm{t}}$ is a false alarm in a Gaussian white noise if the associated patch is not used in the denoising algorithm. In Proposition 3 we choose $a$ in order to control the number of false alarms with high probability.

Proposition 3

Let $\operatorname{NFA}_{\text{max}}\in[0,|T|]$, $T$ given in (10) and let $a\in\mathbb{R}^{\Omega}$ be defined for any $\bm{\mathrm{t}}\in\Omega$ by

[TABLE]

with background model being a Gaussian white noise $W$, i.e. $f=\delta_{0}$ in Definition 4. Let $T$ be defined in (10) and $N_{\omega}(W)\in\{0,\dots,T\}$ the random number of selected patches used to denoise the patch $P_{\omega}(W)$, see Algorithm 2. Then for any $n\in\mathbb{N}\backslash\{0\}$ it holds that

[TABLE]

Proof:

Using the Markov inequality, we have

[TABLE]

$\square$

In this case the null hypothesis $\mathbb{P}_{0}$ is given by a standard Gaussian random field, which is a special case of the Gaussian random field models introduced in Section 3. In the next proposition, using the a contrario framework, we obtain probabilistic guarantees on the distance between the reconstructed patch $\hat{p}(u,\omega)$ and the true patch $P_{\omega}(u_{0})$.

Proposition 4

Let $U=u_{0}+\sigma W$, where $W$ is a standard Gaussian white noise over $\Omega$, $u_{0}\in\mathbb{R}^{\Omega}$ and $\sigma>0$. Let $\bm{\mathrm{x}}\in\Omega$ and $\omega=\bm{\mathrm{x}}+\omega_{0}$ be a fixed patch and let $\operatorname{NFA}_{\text{max}}\in[0,|T|]$. We introduce the random set $\hat{T}=\{\bm{\mathrm{t}}\in T,\ \mathcal{AS}(U,\bm{\mathrm{t}},\omega)\leq\sigma^{2}a(\bm{\mathrm{t}})\}$ (the selected offsets) with $a(\bm{\mathrm{t}})=\mathsf{AP}^{-1}\left(\bm{\mathrm{t}},\omega,1-\operatorname{NFA}_{\text{max}}/|T|\right)$ as in Proposition 3 and $T$ defined in (10). Let $a_{T}=\max_{\bm{\mathrm{t}}\in T}a(\bm{\mathrm{t}})$. Then for any $a_{W}>0$, setting $\varepsilon_{W}=1-\mathbb{P}\left[\|P_{\omega}(W)\|_{2}^{2}\leq a_{W}\ |\ \hat{T}\right]$, we have

[TABLE]

Proof:

We have for any $\bm{\mathrm{t}}\in\hat{T}$

[TABLE]

This gives the following event inclusion for any $\bm{\mathrm{t}}\in\hat{T}$,

[TABLE]

We also have that by definition of $\varepsilon_{W}$

[TABLE]

$\square$

In our applications we use Algorithm 2 with $a(\bm{\mathrm{t}})=\mathsf{AP}^{-1}\left(\bm{\mathrm{t}},\omega,1-\operatorname{NFA}_{\text{max}}/|T|\right)$. Therefore we need to compute $a(\bm{\mathrm{t}})=\mathsf{AP}^{-1}\left(\bm{\mathrm{t}},\omega,1-\operatorname{NFA}_{\text{max}}/|T|\right)$ with a Gaussian white noise background model. We recall that in Section 3.2, using Proposition 2, we give a method to compute this quantity in general Gaussian settings. In the case of a Gaussian white noise, the next proposition shows that the eigenvalues can be computed without approximation.

Proposition 5

Let $\bm{\mathrm{t}}=(t_{x},t_{y})\in\mathbb{Z}^{2}\backslash\{0\}$, $C_{\bm{\mathrm{t}}}$ as in Proposition 2 with $f=\delta_{0}$ and $\omega=\llbracket 0,p-1\rrbracket^{2}$, with $p\in\mathbb{N}$. We have, expressing $C_{\bm{\mathrm{t}}}$ in the basis corresponding to the raster scan order on the $x$-axis

[TABLE]

where $D_{j}$ is a zero matrix with ones on the $j$-th diagonal. The eigenvalues of $C_{\bm{\mathrm{t}}}$ are given by $\lambda_{m,k}=4\sin^{2}\left(\frac{k\pi}{2m}\right)$ with multiplicity $r_{m,k}$ where $m\in\llbracket 2,q+1\rrbracket$, $k\in\llbracket 1,m-1\rrbracket$ and $q=\lceil\frac{p}{|t_{x}|\vee|t_{y}|}\rceil$. For any $m\in\llbracket 2,q+1\rrbracket$, $k\in\llbracket 1,m-1\rrbracket$ it holds

1. (a)

for any $k^{\prime}\in\llbracket 1,m-1\rrbracket$, $r_{m,k}=r_{m,k^{\prime}}\;;$ 2. (b)

$r_{m,k}=2|t_{x}||t_{y}|\ \text{if}\ 2\leq m<q\;;$** 3. (c)

$r_{m,k}=r_{x}r_{y}\ \text{if}\ m=q+1\;;$** 4. (d)

$\sum_{m=2}^{q+1}\sum_{k=1}^{m-1}r_{m,k}=p^{2}\;,$**

with $r_{x}=\left(\lceil\frac{p}{|t_{x}|}\rceil-q\right)|t_{x}|+|t_{x}|-p_{x}$, where $p_{x}=|t_{x}|\lceil\frac{p}{|t_{x}|}\rceil-p$. We define $r_{y}$ in the same manner. A similar proposition holds if $t_{y}\neq 0$.

Proof:

The proof is postponed to Appendix A. $\square$

This property allows us to compute exactly the eigenvalues appearing in Proposition 2. In Figure 3 we illustrate that $a(\bm{\mathrm{t}})$ for fixed patch size ($8\times 8$) and patch search window ($21\times 21$). Thus in our implementation we suppose that $a(\bm{\mathrm{t}})=\mathsf{AP}^{-1}\left(\bm{\mathrm{t}},\omega,1-\operatorname{NFA}_{\text{max}}/|T|\right)$ is constant and set its value to the mean of $a(\bm{\mathrm{t}})$ over $\bm{\mathrm{t}}\in T$.

4.2 Some experimental results

In the following paragraph we present and comment some results of our threshold NL-means algorithm, see Algorithm 2. We recall that we use $a(\bm{\mathrm{t}})=\sum_{\bm{\mathrm{t}}\in T}\mathsf{AP}^{-1}\left(\bm{\mathrm{t}},\omega,1-\operatorname{NFA}_{\text{max}}/|T|\right)/|T|$. In Figure 4 we present a first comparison with the NL-means algorithm. Perceptual results as well as Peak Signal to Noise Ratio ($\operatorname{PSNR}$) measurements 111$\operatorname{PNSR}(u,v)=10\log_{10}\left(\frac{\max_{\Omega}u^{2}}{\|u-v\|_{2}^{2}}\right)\;.$ are commented. We also present the running time of the original NL-means algorithm and ours. The experiments were conducted with the following computer specifications: 16G RAM, 4 Intel Core i7-7500U CPU (2.70GHz). Results on other images than Barbara are displayed in Figure 5.

If the threshold $a(\bm{\mathrm{t}})$ is high, i.e. $\operatorname{NFA}_{\text{max}}\ll|T|$ then almost no patch is rejected, which means that almost all patches are used in the denoising process. In consequence, the output denoised image is very smooth. This smoothness is a correct guess for constant patches. However, this proposition does not hold when the region contains details. Indeed, in this case details are lost due to the averaging process. By setting a conservative threshold, e.g. $\operatorname{NFA}_{\text{max}}/|T|\approx 0.1$, for example, we reject all the patches for which the structure does not strongly match the one of the input patch, see Figure 6. This conservative property of the algorithm ensures that we can control the loss of information in the denoised image, see Proposition 4. However, if no patch, other than the input patch itself, is detected as similar we highly overfit the original noise. Many algorithms such as BM3D, see [12], solve this problem by treating this case as an exception, applying a specific denoising method in this situation. We show the differences between our version of NL-means and BM3D in Figure 7 .

In Figure 8, we show that Algorithm 2 performs better than the original NL-means algorithm. By setting $\operatorname{NFA}_{\text{max}}/|T|=0.01$ we obtain that the $\operatorname{PSNR}$ of the denoised image is better than the one of NL-means for nearly every value of $h$.

Let us emphasize that our goal is not to provide a new state-of-the-art denoising algorithm. Indeed we never obtain better denoising results than the BM3D algorithm. However, our algorithm slightly improves the original NL-means algorithm. It shows that statistical testing can be efficiently used to measure the similarity between patches and therefore provides a robust way to perform the weighted average in this algorithm.

5 Periodicity analysis

5.1 Existing algorithms

In the following sections we use our patch similarity detection algorithm, see Algorithm 1, to analyze images exhibiting periodicity features. Let $\Omega\subset\mathbb{Z}^{2}$ be a finite domain and $\omega\subset\Omega$ a finite patch domain.

Periodicity detection is a long-standing problem in texture analysis [71]. First algorithms used the quantization of images, relying on co-occurrence matrices and statistical tools like $\chi_{2}$ tests or $\kappa$ tests. Global methods extract peaks in the frequency domain (Fourier spectrum) [49] or in the spatial domain (autocorrelation). In [32] the notion of inertia is introduced. It is defined for any $\bm{\mathrm{t}}\in\Omega$ by $\mathcal{I}(\bm{\mathrm{t}})=\sum_{(i,j)\in\llbracket 0,N_{g}\rrbracket^{2}}(i-j)^{2}\left(\sum_{\bm{\mathrm{z}}\in\Omega}\mathbb{1}_{\dot{u}(\bm{\mathrm{z}})=i}\mathbb{1}_{\dot{u}(\bm{\mathrm{z+t}})=j}\right)$, where $u$ is a quantized image on $N_{g}+1$ gray levels. In [8], the authors show that the local minima of the inertia measurement can be used to assess periodicity. Similarly, we introduce the $\omega$-inertia for any $\bm{\mathrm{t}}\in\Omega$ by $\mathcal{I}_{\omega}(\bm{\mathrm{t}})=\sum_{(i,j)\in\llbracket 0,N_{g}\rrbracket^{2}}(i-j)^{2}\left(\sum_{\bm{\mathrm{z}}\in\omega}\mathbb{1}_{\dot{u}(\bm{\mathrm{z}})=i}\mathbb{1}_{\dot{u}(\bm{\mathrm{z+t}})=j}\right)$. The following proposition extends to a local framework results from [52].

Proposition 6

Let $u\in\mathbb{R}^{\Omega}$. Suppose that $u$ is quantized, i.e. there exists $N_{g}\in\mathbb{N}$ such that for any $\bm{\mathrm{x}}\in\Omega$, $u(\bm{\mathrm{x}})\in\llbracket 0,N_{g}\rrbracket$. We have $\mathcal{I}_{\omega}(\bm{\mathrm{t}})=\mathcal{AS}(u,\bm{\mathrm{t}},\omega)$.

Proof:

For any $\bm{\mathrm{t}}\in\Omega$ we have

[TABLE]

$\square$

If $\omega=\Omega$ then the $\omega$-inertia statistics is exactly the inertia introduced in [32] and the result is due to [52].

5.2 Algorithm and properties

Lattice detection is closely related to periodicity analysis, since identifying a lattice is similar to extracting periodic or pseudo-periodic structures up to deformations and approximations. A state-of-the-art algorithm proposed in [53] uses a recursive framework which consists in 1) a lattice model proposal based on detectors such as Kanade-Lucas-Tomasi ($\operatorname{KLT}$) feature trackers [48], 2) spatial tracking using inference in a probabilistic graphical model, 3) spatial warping correcting the lattice deformations in the original image. In this section we propose a new algorithm for lattice detection. The lattice proposal step 1) is replaced by an Euclidean auto-similarity matching detection (see Section 3.2 and Algorithm 1) where the patch domain $\omega$ is fixed. Using these detections we build a graph with a few nodes (usually approximately $20$ nodes for a $256\times 256$ image). We use the same notation for the detection mapping $\bm{\mathrm{t}}\mapsto\mathbb{1}_{\mathcal{AS}_{i}(u,\bm{\mathrm{t}},\omega)\leq a(\bm{\mathrm{t}})}$ i.e. the $D_{map}$ output of Algorithm 1, which is a binary function over the offsets, and the set of detected offsets. We recall that two pixel coordinates $\bm{\mathrm{x}}$ and $\bm{\mathrm{y}}$ are said to be $8$-connected if $\bm{\mathrm{x}}=\bm{\mathrm{y}}+(\delta_{x},\delta_{y})$ with $\delta_{x},\delta_{y}\in\{-1,0,1\}$. The graph $\mathscr{G}=(V,E)$ is then built as follows:

• $\blacktriangleright$

Vertices: for each 8-connected component, $\mathscr{C}_{k}$ in $D_{map}$ we note $\bm{\mathrm{v}}_{k}$ one position for which the minimum of $\mathcal{AS}(u,\bm{\mathrm{t}},\omega)$ over $\mathscr{C}_{k}$ is achieved. The set of vertices $V$ is defined as $V=\left(\bm{\mathrm{v}}_{k}\right)_{k\in\llbracket 1,N_{\mathscr{C}}\rrbracket}$ where $N_{\mathscr{C}}$ is the number of 8-connected components in $D_{map}$ ;

• $\blacktriangleright$

Edges: each vertex is linked with its four nearest neighbors in the sense of the Euclidean distance, defining four unoriented edges.

Refering to the three steps of [53] we present our model to replace step 2) (i.e. the inference in a probabilistic graphical model) and introduce the approximated lattice hypothesis defined on a graph.

Definition 5** (Approximated lattice hypothesis)**

Let $\mathscr{G}=(V,E)$ be a random graph with $V\subset\mathbb{R}^{2}$. We say that $\mathscr{G}$ follows the approximated lattice hypothesis if there exists a basis $B=(b_{1},b_{2})$ of $\mathbb{R}^{2}$ and, for each edge $\bm{\mathrm{e}}\in E$, a couple of integers $(m_{\bm{\mathrm{e}}},n_{\bm{\mathrm{e}}})\in\mathbb{Z}^{2}$ such that $\bm{\mathrm{e}}$ has the same distribution as $m_{\bm{\mathrm{e}}}b_{1}+n_{\bm{\mathrm{e}}}b_{2}+\sigma Z_{\bm{\mathrm{e}}}$, with $(Z_{\bm{\mathrm{e}}})_{\bm{\mathrm{e}}\in E}$ independent standard Gaussian random variables in $\mathbb{R}^{2}$ and $\sigma>0$. We denote by $M$ the vector of all coefficients $(m_{\bm{\mathrm{e}}},n_{\bm{\mathrm{e}}})_{\bm{\mathrm{e}}\in E}\in\mathbb{Z}^{2|E|}$.

Our goal is to perform inference in the statistical model defined by the following log-posterior

[TABLE]

where $r(B,M)=\delta_{B}\|B\|_{2}^{2}+\delta_{M}\|M\|_{2}^{2}$ with $\delta_{B},\delta_{M}\geq 0$. A discussion on the dependence of the model on the hyperparameters $(\delta_{B},\delta_{M})$ is conducted in Figure 9. Finding the $\operatorname{MLE}$ of this full log-posterior is a non-convex, integer problem. However performing the minimization alternatively on $B$ and $M$ is easier since at each step we only have a quadratic function to minimize. Minimizing a positive-definite quadratic function over $\mathbb{Z}^{2}$ is equivalent to finding the vector of minimum norm in a lattice. This last formulation is known as the Shortest Vector Problem ($\operatorname{SVP}$) which is a challenging problem [51] (though it is not known if it is a $\operatorname{NP}$-hard problem). We replace this minimization procedure over a lattice by a minimization over $\mathbb{R}^{2}$ followed by a rounding of this relaxed solution.

For any $\sigma>0$ we denote by $\mathscr{L}_{n}(\sigma)=\mathscr{L}(B_{n},M_{n},\sigma^{2}|E)$, with $n\in\mathbb{N}$, the log-posterior sequence.

Proposition 7** (Alternate minimization update rule)**

In Algorithm 3, we get for any $n\in\mathbb{N}$

[TABLE]

with $\otimes$ the tensor product between matrices and

1. (a)

$\Lambda_{B}=\left(\begin{matrix}\|b_{1}\|^{2}+\delta_{B}&\langle b_{1},b_{2}\rangle\\ \langle b_{1},b_{2}\rangle&\|b_{2}\|^{2}+\delta_{B}\end{matrix}\right)\;,\qquad\Lambda_{M}=\left(\begin{matrix}\|M_{1}\|^{2}+\delta_{M}&\langle M_{1},M_{2}\rangle\\ \langle M_{1},M_{2}\rangle&\|M_{2}\|^{2}+\delta_{M}\end{matrix}\right)\;;$** 2. (b)

$E_{B}=\left(\begin{matrix}(\langle\bm{\mathrm{e}},b_{1}\rangle)_{\bm{\mathrm{e}}\in E}\\ (\langle\bm{\mathrm{e}},b_{2}\rangle)_{\bm{\mathrm{e}}\in E}\end{matrix}\right)\;,\qquad E_{M}=\left(\begin{matrix}\underset{\bm{\mathrm{e}}\in E}{\overset{}{\sum}}{m_{\bm{\mathrm{e}}}\bm{\mathrm{e}}}\\ \underset{\bm{\mathrm{e}}\in E}{\overset{}{\sum}}{n_{\bm{\mathrm{e}}}\bm{\mathrm{e}}}\end{matrix}\right)\;.$**

Proof:

The proof is postponed to Appendix B. $\square$

Note that if $B$ is orthogonal, i.e. $\langle b_{1},b_{2}\rangle=0$ then $\Lambda_{B}$ is diagonal and the proposed method is the exact solution to the minimization problem over $\mathbb{Z}^{2}$.

Theorem 1** (Convergence in finite time)**

For any $\sigma>0$, $(\mathscr{L}_{n}(\sigma))_{n\in\mathbb{N}}$ is a non-decreasing sequence. In addition, $\left(B_{n}\right)_{n\in\mathbb{N}}$ and $\left(M_{n}\right)_{n\in\mathbb{N}}$ converge in a finite number of iterations.

Proof:

$(\mathscr{L}_{n}(\sigma))_{n\in\mathbb{N}}$ is non-decreasing since for any $n\in\mathbb{N}$, $\mathscr{L}_{n}(\sigma)\leq\mathscr{L}(B_{n},M_{n+1},\sigma^{2}|E)\leq\mathscr{L}_{n+1}(\sigma)$. Let us show that the sequences $(M_{n})_{n\in\mathbb{N}}$ and $(B_{n})_{n\in\mathbb{N}}$ are bounded. Because $(\mathscr{L}_{n}(\sigma))_{n\in\mathbb{N}}$ is non-decreasing, the sequence $\left(q(B_{n},M_{n}|E)\right)_{n\in\mathbb{N}}$ is non-increasing. We obtain that

[TABLE]

The sequence $\left(M_{n}\right)_{n\in\mathbb{N}}$ is bounded thus we can extract a converging subsequence. Since $\left(M_{n}\right)_{n\in\mathbb{N}}$ takes value in $\mathbb{Z}^{2|E|}$, this subsequence is stationary with value $M$. Let $n_{0}\in\mathbb{N}$ be the first time we hit value $M$. Let $n\in\mathbb{N}$, with $n\geq n_{0}+1$, there exists $n_{1}\in\mathbb{N}$, with $n_{1}\geq n$ such that $M_{n_{1}}=M_{n_{0}}$ thus

[TABLE]

Hence for every $n\geq n_{0}+1$, $\mathscr{L}_{n}(\sigma)=\mathscr{L}(B_{n},M_{n},\sigma^{2}|E)=\tilde{\mathscr{L}}(\sigma)$. Suppose there exists $n\geq n_{0}+1$ such that $M_{n}\neq M_{n+1}$ this means that $\mathscr{L}(B_{n},M_{n+1},\sigma^{2}|E)>\mathscr{L}_{n}(\sigma)$ (because of lines 6 and 7 of Algorithm 3) which is absurd. Thus $\left(M_{n}\right)_{n\in\mathbb{N}}$ is stationary and so is $\left(B_{n}\right)_{n\in\mathbb{N}}$. $\square$

In Algorithm 3 $M_{0}$ is initialized with zero and $B_{0}$ is defined as an orthonormal (up to a dilatation factor) direct basis where the first vector is given by an edge with median norm in $E$.

5.3 Experimental results

Combining the results of Section 5.2 and Section 3.2 we obtain an algorithm to extract lattices in images, see Figure 10. In what follows we perform lattice detection using Algorithm 1 in order to extract auto-similarity given a patch in an original image $u$, which implies that the patch domain $\omega$ is set by the user. Recall that in Algorithm 1, the eigenvalues of the covariance matrix in Proposition 2 are approximated, and that the cumulative distribution function of the quadratic form in Gaussian random variables is computed via the Wood F method [66]. Lattice detection is performed using Algorithm 3 with parameters $\delta_{M}=10$ and $\delta_{B}=10^{-2}$.

5.3.1 Escher paving

In this section we study art images, Escher pavings, with strongly periodic structure. We investigate the following parameters of our lattice detection algorithm:

1. (a)

background microtexture model $\mathbb{P}_{0}$, 2. (b)

$\operatorname{NFA}_{\text{max}}$ parameter in Algorithm 1, 3. (c)

patch domain $\omega$.

Microtexture model

We confirm that the choice of the microtexture model will influence the detected geometrical structures. The more structured is the background noise model the less we obtain detections. This situation is considered in Figure 11.

$\operatorname{NFA}_{\text{max}}$ parameter

Using a more adapted microtexture model as background model we gain robustness compared to other less structured models such as a Gaussian white noise. However, $\operatorname{NFA}_{\text{max}}$ must be set carefully otherwise two situations can occur:

1. (a)

if $\operatorname{NFA}_{\text{max}}$ is too high, too many detections can be obtained (true perceptual detections are not differentiated from false positives) ; 2. (b)

if $\operatorname{NFA}_{\text{max}}$ is too low, we fail to identify important perceptual structures in the image.

We observe that a general good practice is to set $\operatorname{NFA}_{\text{max}}$ equal to $10$, see Figure 12. However, if the input patch is corrupted one may increase this parameter up to $10^{2}$ or $10^{3}$, see Figure 17 and Figure 18.

Patch position

Patch position and size are crucial in our detection model, since we rely on local properties of the image. As shown in Figure 13 these parameters should be carefully selected by the user. However, for particular applications such as lattice extraction for crystallographic purposes, there exist procedures to extract primitive cells [50].

5.3.2 Crystallography images

Defect localization, noise reduction, correction of crystalline structures in images are central tasks in crystallography. Usually, they require the knowledge of the geometry of a perfect underlying crystal. In our experiments we manually identify the geometry of the periodic crystal, which allows for multiple structures in one image, provided a user input of the primitive cell in a lattice. This primitive cell extraction could be automated [50]. In Figure 14, we present an example of multiple geometry extraction. Statistics like angle and period can be retrieved using the estimated basis vectors. This image contains two lattices and the locality of our measurements allows for the detection of both structures. Using windowed Fourier transform could be efficient to obtain local measurements on the periodicity of these images since the information is highly frequential. However in order to obtain the same detection map as Algorithm 1 one must carefully set the threshold parameter, $\operatorname{NFA}_{\text{max}}$. This situation is illustrated in Figure 15.

Finally we assess the precision of our measurements by comparing our results with a model used in crystallography, see Figure 16. We indeed retrieve one of the possible bases used to describe these lattices. However, the symmetry constraints are not present in the identified basis. To obtain another basis, one must relax the regularization parameters. A more natural way to obtain the desired primitive cell would be to introduce symmetry constraints in the graphical model formulation in (14).

5.3.3 Natural images

Identifying lattices in natural images is a more challenging task since we have to deal with image artifacts. In this section we investigate the effect on the detection of the background clutter in natural images, see Figure 17, and the effect of the camera position, see Figure 18.

Preprocessing

Due to the occlusions occurring in natural images, if a lattice is superposed over a real photograph, carefully selecting structural elements might not be enough in order to retrieve the periodicity. Indeed, if we observe a repetition of the lattice pattern, the background does not necessarily contain any repetition and thus makes the detection more complicated. In order to avoid such a problem we propose to introduce a preprocessing step in our algorithm. This preprocessing step will be encoded in a linear filter $h$. Suppose $U$ is a sample from a Gaussian model with function $f$ then $h*U$ is a sample from a Gaussian model with function $h*f$. Thus all the properties derived earlier remain valid with this linear operation. In Figure 17, we set $h$ to be a Laplacian operator 222We use a discrete Laplacian operator $\Delta$ such that for any $\bm{\mathrm{x}}=(x_{1},x_{2})$, we get that $\Delta(u)(x_{1},x_{2})=\left(u(x_{1}+1,x_{2})+u(x_{1}-1,x_{2})+u(x_{1},x_{2}+1)+u(x_{1},x_{2}-1)-4u(\bm{\mathrm{x}})\right)/4$, where boundaries are handled periodically.. This operation allows us to avoid contrast problems.

Homography

In the previous experiments we suppose that the lattice structure was in front of the camera. In many cases this assumption is not true and there exists an homography that matches the deformed lattice in the image to a true lattice. Our algorithm makes the assumption that the lattice is viewed in a frontal way and fails otherwise. However, locally, this assumption is true and we can observe partial match of the lattices in Figure 18.

5.4 Texture ranking

We conclude these experiments by showing that this simple graphical model can be used to perform ranking among texture images, sorting them according to their degree of periodicity. We say that an image has high periodicity degree if a lattice structure can be well fitted to the image. We introduce a criterion for evaluating the relevance of the lattice hypothesis. Let $u$ be an image over $\Omega$, let $\omega\subset\Omega$ be a patch domain and $a$ be as in Proposition 1 with $\operatorname{NFA}_{\text{max}}$ set by the user.

Definition 6** (Periodicity criterion)**

Let $\{\bm{\mathrm{t}}\in\Omega,\ \mathcal{AS}(u,\bm{\mathrm{t}},\omega)\leq a(\bm{\mathrm{t}})\}$ be the set of detected offsets and $N_{{\mathscr{C}}}$ its number of connected components as defined in Section 5.2. Let also $(\widehat{B},\widehat{M},\widehat{\sigma})$ be the estimated parameters using Algorithm 3. We define the following periodicity criterion $c_{per}$ as

[TABLE]

where $\widehat{B}=(\hat{b}_{1},\hat{b}_{2})$.

The criterion $c_{per}$ simply computes the ratio between the error area of Algorithm 3, i.e. the error made when considering the approximated lattice hypothesis, see Definition 5, and the area of the parallelogram defined by the output basis vectors. If we have enough detections this quantity is supposed to be small when the approximated lattice hypothesis holds and large when it does not. Nonetheless, we introduce a dependency in the number of detections. Indeed, even if no lattice is perceived, the hypothesis in Definition 5 may still hold if the number of detected offsets is small.

In the experiment presented in Figure 19 we sort 25 texture images based on the $c_{per}$ criterion. Images are of size ${256\times 256}$. Since the identified graph highly depends on the patch position and the patch size, for each image we uniformly sample 150 patch positions and set the patch size to ${20\times 20}$. For each set of parameters we find a lattice using Algorithm 1 and Algorithm 3 with parameters \text{\operatorname{NFA}_{\text{max}}}=1, $\delta_{M}=10$, $\delta_{B}=10^{-2}$ and $N_{it}=10$. A statistical study of our ranking is presented in Figure 20. Note that, from a perceptual point of view, from (a) to (n) all textures are periodic except for (f), (j) and (k) which are examples for which our algorithm fails. However, from (o) to (y), no texture is periodic.

6 Conclusion

In this paper we introduce a statistical model, the a contrario framework, to analyze spatial redundancy in images. We propose a general algorithm for detecting redundancy in natural images. It relies on Gaussian random fields as background models and takes advantage of the links between the $\ell^{2}$ norm and Gaussian densities. The a contrario formulation provides us with a statistically sound way of thresholding distances in order to assess similarity between patches. In this rationale we replace the task of manually setting thresholds by the selection of a Number of False Alarms.

We illustrate our contribution with three examples in various domains of image processing. Introducing a simple modification of the NL-means algorithm we show that similarity detection (in this case, dissimilarity detection) in a theoretical a contrario framework can easily be embedded in any image denoising pipeline. For instance, the threshold we introduced could be integrated into the Non-Local Bayes algorithm [41] in order to estimate mean and covariance matrices with probabilistic guarantees. The generality of our model allows for several extensions for non-Gaussian noises [16] or to take into account the geometry of the patch space [34, 63].

Turning to periodicity detection we propose a novel graphical model using the output of Algorithm 1 in order to extract lattices from images. In this model, lattice extraction is formulated as the maximization of some log-likelihood defined on a graph. We prove the finite-time convergence of Algorithm 3. We provide image experiments illustrating the role of the hyperparameters in our study and we assess the importance of selecting adaptive Gaussian random fields as background models. We remark that the expected number of false alarms parameter is linked to the choice of the input patch and give a range of possible values for $\operatorname{NFA}_{\text{max}}$ settings. We also illustrate its possible application in crystallography as it correctly identifies underlying lattices in alloys. This rationale could be used to identify symmetry groups (wallpaper groups) in alloys, following the work of [45]. Finally our method is tested on natural images where some of its limits such as perspective defect or sensitivity to occlusion phenoma are identified. It must be noted that our method could easily be extended to color images by considering $\mathbb{R}^{3}$-valued instead of real-valued images.

Our last application consists in giving a quantitative criterion for periodicity texture ranking. This criterion is based on the parameters estimated in Algorithm 3. Since we set our background models to be Gaussian random fields and remarking that these are good microtexture approximations we wish to explore the possibility to embed our a contrario framework in texture analysis and texture synthesis algorithms. For instance an a contrario methodology could be incorporated in the algorithm proposed by Raad et al. in [56]. Another potential direction is to look at the behavior of the introduced dissimilarity functions for more general random fields in order to handle more complex and structured situations such as parametric texture synthesis.

7 Acknowledgements

The authors would like to thank Denis Gratias for the crystallography images, Jérémy Anger for some of natural images, Axel Davy who provided an OpenCL implementation of the NL-means algorithm and Thibaud Ehret for its insights and comments on denoising algorithms.

Appendix A Eigenvalues

Proof:

[Proof of Proposition 5]

We fix $\bm{\mathrm{t}}\neq 0$ with $\|\bm{\mathrm{t}}\|_{\infty}<p$ and denote $C=C_{\bm{\mathrm{t}}}$. Without loss of generality we consider that $t_{x}>0$ and $t_{y}>0$. We consider $X$ an eigenvector of $C$. Let $\Omega_{0}=\left(\Omega-\bm{\mathrm{t}}\right)\cap\Omega^{c}$ and the function $J:\ \Omega_{0}\to\llbracket 2,+\infty\llbracket$ such that for any $\bm{\mathrm{x}}_{0}\in\Omega_{0}$

[TABLE]

It is clear that $I=\{(k,m),k\in\llbracket 1,m-1\rrbracket,\ m\in J(\Omega_{0})\}$ is in bijection with $\Omega$. Let $\bm{\mathrm{x}}_{0}\in\Omega_{0}$, $m=J(\bm{\mathrm{x}}_{0})$ and $k\in\llbracket 1,m-1\rrbracket$. We define $X_{\bm{\mathrm{x}}_{0},k}$ over $\mathbb{Z}^{2}$ such that

[TABLE]

Using that $\sin(a+b)+\sin(a-b)=2\sin(a)\cos(b)$, we have for any $\bm{\mathrm{x}}\in\mathbb{Z}^{2}$

[TABLE]

This implies that for any $\bm{\mathrm{x}}\in\mathbb{Z}^{2}$

[TABLE]

Thus the one-dimensional vector (given by the raster-scan order on the $x$-axis) of the restriction of $X_{\bm{\mathrm{x}}_{0},k}$ is an eigenvector of $C$ associated with eigenvalue $4\sin^{2}\left(\frac{k\pi}{m}\right)$.

Using that $I$ is in bijection with $\Omega$ we get that the number of vectors $(X_{\bm{\mathrm{x}}_{0},k})$ is $|\Omega|$. We show that this family of vectors is linearly-independent. Let $\Lambda_{\bm{\mathrm{x}}_{0},k}\in\mathbb{R}$ such that

[TABLE]

Since $X_{\bm{\mathrm{x}}_{0},k}$ and $X_{\bm{\mathrm{y}}_{0},k^{\prime}}$ have different support if and only if $\bm{\mathrm{x}}_{0}\neq\bm{\mathrm{y}}_{0}$ we get that for any $\bm{\mathrm{x}}_{0}\in\Omega_{0}$, $\sum_{k=1}^{J(\bm{\mathrm{x}}_{0})-1}\Lambda_{\bm{\mathrm{x}}_{0},k}X_{\bm{\mathrm{x}}_{0},k}=0$. This gives that $(\Lambda_{\bm{\mathrm{x}}_{0},k})_{k\in\llbracket 1,J(\bm{\mathrm{x}}_{0})-1\rrbracket}$ is in the kernel of the matrix $\left(\sin(\ell k\pi/(J(\bm{\mathrm{x}}_{0})-1))\right)_{1\leq j,\ell\leq J(\bm{\mathrm{x}}_{0})-1}$. Since the sinus discrete transform is invertible we obtain that for any $\bm{\mathrm{x}}_{0}\in\Omega_{0}$ and $k\in\llbracket 1,J(\bm{\mathrm{x}}_{0})-1\rrbracket$, $\Lambda_{\bm{\mathrm{x}}_{0},k}=0$. Thus the family $X_{\bm{\mathrm{x}}_{0},k}$ is a basis of eigenvectors.

We aim at computing the cardinality of $K_{k,m}=\{X_{\bm{\mathrm{x}}_{0},k},J(\bm{\mathrm{x}}_{0})=m\}$. By definition, in Proposition 5, $r_{k,m}=|K_{k,m}|$. First note that $|K_{k^{\prime},m}|=|K_{k,m}|$. We give the following decomposition $\Omega_{0}=\Omega_{x}\cup\Omega_{y}\cup\Omega_{x,y}$ with

[TABLE]

Note that for all $\bm{\mathrm{x}}_{0}\in\Omega_{0}$ we have that $\bm{\mathrm{x}}_{0}+(q+1)\bm{\mathrm{t}}\notin\Omega$, with $q=\lceil\frac{p}{|t_{x}|\vee|t_{y}|}\rceil$. Thus $J(\Omega_{0})\subset\llbracket 2,q+1\rrbracket$. Let $m\in\llbracket 2,q-1\rrbracket$. The cardinality of $K_{k,m}$ is the cardinality of $J^{-1}(m)$. Let $\bm{\mathrm{x}}_{0}\in\Omega_{x}$ we have

[TABLE]

Since $\bm{\mathrm{x}}_{0}\in\Omega_{x}$ we have $i_{0}+mtx\leq p-1$, hence

[TABLE]

Thus $|\Omega_{x}\cap J^{-1}(m)|=t_{x}t_{y}$. Similarly we get that $|\Omega_{y}\cap J^{-1}(m)|=t_{x}t_{y}$ and $\Omega_{x,y}\cap J^{-1}(m)=\emptyset$. Thus, $|K_{k,m}|=2t_{x}t_{y}$.

We have computed $|K_{k,m}|$ for every $m\in\llbracket 2,q-1\rrbracket$. In order to complete our study it only remains to compute $|K_{k,q+1}|$, since $|K_{k,q}|$ can be deduced from the summability condition and from $|K_{k,m}|=|K_{k^{\prime},m}|$. We only compute $|K_{k,q+1}|$. We remark that $\Omega_{x}\cap J^{-1}(q+1)=\Omega_{y}\cap J^{-1}(q+1)=\emptyset$. Let $\bm{\mathrm{x}}_{0}\in\Omega_{x,y}$ then $\bm{\mathrm{x}}_{0}=-\bm{\mathrm{t}}+(x,y)$ with $x\in\llbracket 0,t_{x}-1\rrbracket$ and $y\in\llbracket 0,t_{y}-1\rrbracket$. We obtain the following equivalence

[TABLE]

Since $qt_{x}\geq p$ or $qt_{y}\geq p$ we obtain that the first condition is always satisfied. Thus we get

[TABLE]

Using that $p-1-(q-1)t_{x}=\left(\lceil\frac{p}{t_{x}}\rceil-q\right)t_{x}+t_{x}-1-p_{x}$, we conclude the proof. $\square$

Appendix B Update rules

We derive the proof of Proposition 7.

Proof:

Computing the minimum of $q(B,M|E)$ for fixed $B\in\mathbb{R}^{4}$, respectively fixed $M\in\mathbb{R}^{2|E|}$, gives the update rule for $M$, respectively for $B$. We obtain that

[TABLE]

where $\alpha(M)$ depends only on $M$. Similar derivation goes for $B$ and we obtain the proposed update rules. $\square$