Redundancy in Gaussian random fields

In this paper, we introduce a notion of spatial redundancy in Gaussian random fields. This study is motivated by applications of the a contrario method in image processing. We define similarity functions on local windows in random fields over discrete or continuous domains. We derive explicit Gaussian asymptotics for the distribution of similarity functions when computed on Gaussian random fields. Moreover, for the special case of the squared L2 norm, we give non-asymptotic expressions in both discrete and continuous periodic settings. Finally, we present fast and accurate approximations of these non-asymptotic expressions using moment methods and matrix projections.

random fields [26,40,55,57] give good visual results for input textures with no, or few, spatial redundancy, they fail when it comes to sampling structured textures (brick walls, fabric with repeated patterns, etc.). In this case, more elaborated models are needed [20,27,43]. In this work, we derive explicit probability distribution functions for the random variables associated with the output of similarity functions computed on local windows of random fields. The knowledge of such functions allows us to conduct rigorous statistical testing on the spatial redundancy in natural images.
In order to compute these explicit distributions we will consider specific random fields over specific topological spaces. First, the random fields will be defined either over R 2 (or T 2 , where T 2 is the 2-dimensional torus, when considering periodicity assumptions on the field), or over Z 2 (or pZ{pM Zqq 2 , with M P N when considering periodicity assumptions on the field). Each of these spaces is embedded with its classical topology. The first case is the continuous setting, whereas the second one is the discrete setting. In image processing, the most common framework is the finite discrete setting. The discrete setting (Z 2 ) can be used to define asymptotic properties when the size of images grows or when their resolution increases [8], whereas continuous settings are needed in specific applications where, for instance, rotation invariant models are required [54]. All the considered random fields will be Gaussian. This assumption will allow us to explicitly derive moments of some similarity functions computed on local windows of the random field. Once again, another reason for this restriction comes from image processing. Indeed, given an input image, we can compute its first and second-order statistics. Sampling from the associated Gaussian random field gives examples of images which preserve the covariance structure but lose the global arrangement of the input image. Investigating redundancy of such fields is a first step towards giving a mathematical description of this lack of structure.
Finding measurements which correspond to the ones of our visual system is a long-standing problem in image processing. It was considered in the early days of texture synthesis and analyzed by Julesz [36,37,58] who formulated the conjecture that textures with similar first-order statistics (first conjecture) or that textures with similar first and second-order statistics (second conjecture) could not be discriminated by the human eye. Even if both conjectures were disproved [22], the work of Gatys et al. [27] suggests that second-order statistics of image features are enough to characterize a broad range of textures. To compute features on images we embed them in a higher dimensional space. This operation can be conducted using linear filtering [47] or convolutional neural networks [27] for instance. Some recent works examine the response of convolutional neural network to elementary geometrical pattern [45], giving insight about the perceptual properties of such a lifting. In the present work, we focus on another embedding given by considering a square neighborhood, called a patch, around each pixel. This embedding, is exploited in many image processing tasks such as inpainting [30], denoising [9,39], texture synthesis [23,24,41,49], etc.
In the special case where the similarity functions are given by the L 2 norm, explicit distributions can be inferred even in the non-asymptotic case. Calculating this distribution exactly is demanding since it requires the knowledge of some covariance matrix eigenvalues as well as an efficient method to compute cumulative distribution functions of quadratic forms of Gaussian random variables. We propose an efficient algorithm to approximate this distribution. In [7], this algorithm is applied to denoising and periodicity detection problems in an a contrario framework.
The paper is organized as follows. We recall basic notions of Gaussian random fields in general settings in Section 2.1. Similarity functions to be evaluated on these random fields, as well as their statistical properties, are described in Section 2.2. We give the asymptotic properties of these similarity functions in Gaussian random fields in the discrete setting in Section 3.1 and in the continuous setting in Section 3.2. It is shown in Section 3.3 that the Gaussian asymptotic approximation is valid only for large patches. In order to overcome this problem we consider an explicit formulation of the probability distribution function for a particular similarity function: the square L 2 norm. The computations are conducted in the finite discrete case in Section 4.1. We also derive an efficient algorithm to compute these probability distribution functions. Similar non-asymptotic expressions are given in the continuous case in Section 4.2. Technical proofs and additional results on multidimensional central limit theorems are presented in the Appendices. Lemma 2.3 (Sample path continuity). Assume (A2.1) or (A2.2). In addition, assume that for any x P Ω, mpxq " 0. Then there exists a modification of U , i.e. a random field r U such that for any x P Ω, PrU pxq " r U pxqs " 1, and for any a P A, r U paq is continuous over Ω.
In the rest of the paper we always replace U by its continuous modification r U . Note that in the discrete case all random fields are continuous with respect to the discrete topology.
In Sections 3 and 4, we will suppose that U is a stationary Gaussian random field with zero mean. Asymptotic theorems derived in the next section remain true in broader frameworks, however restricting ourselves to stationary Gaussian random fields allows for explicit computations of asymptotic quantities in order to numerically assess the rate of convergence.

Similarity functions
In order to evaluate redundancy in random fields, we first need to derive a criterion for comparing random fields. We introduce similarity functions which take rectangular restrictions of random fields as inputs.
When comparing local windows of random fields (patches), two cases can occur. We can compare a patch with a patch extracted from the same image. We call this situation internal matching. Applications can be found in denoising [9] or inpainting [17] where the information of the image itself is used to perform the image processing task. On the other hand, we can compare a patch with a patch extracted from another image. We call this situation template matching. An application can be found in the non-parametric exemplar-based texture synthesis algorithm proposed by Efros and Leung [23].
The L 2 norm is the usual way to measure the similarity between patches [39] but many other measurements exist, corresponding to different structural properties, see Figure 1.
Depending on the case dx is either the Lebesgue measure or the discrete measure over ω.
The locality of the measurements is ensured by the fact that these functions are defined on patches, i.e. local windows. Following conditions (1) and (3) in [21] one can check that similarity functions (a), (c) and (e) satisfy the following properties -(Symmetry) spP, Qq " spQ, P q; -(Maximal self-similarity) spP, P q ď spP, Qq; -(Equal self-similarities) spP, P q " spQ, Qq.
Note that since s sc , the scalar product similarity, is homogeneous in P , maximal self-similarity and equal selfsimilarity properties are not satisfied. All introduced similarities satisfy the symmetry condition and s 8 satisfies the maximal self-similarity property. In [21], the authors present many other similarity functions all relying on statistical properties such as likelihood ratios, joint likelihood criteria and mutual information kernels. In the present paper we focus only on similarity functions defined directly in the spatial domain.
Definition 2.5 (Auto-similarity and template similarity). Let u and v be two functions defined over a domain Ω Ă R 2 or Z 2 . Let ω Ă Ω be a patch domain. We introduce P ω puq " u| ω , the restriction of u to the patch domain ω. When it is defined we introduce the auto-similarity with patch domain ω and offset t P R 2 or Z 2 Figure 1. Structural properties of similarity functions. In this experiment the image size is 512ˆ512 and the patch size is 20ˆ20. We show the 20 closest patches (red squares) to the upper-left patch (green square) among all patches for different similarity functions. All introduced similarity functions, see Definition 2.4, correctly identify the structure of the patch, i.e. a large clear part with diagonal textures and a dark ray on the right side of the patch, except for s 8 which is too sensitive to outliers. Similarities s 2 , s 1 and s cos have analogous behaviors and identify correct regions. Similarity s sc is too sensitive to contrast and, as it selects a correct patch, it gives too much importance to illumination.
where s i corresponds to s p with p P p0,`8s, s p,p with p P p0,`8q, s sc or s cos . In the same way, when it is defined, we introduce the template similarity with patch ω and offset t by Note that in the finite discrete setting, i.e. Ω " pZ{pM Zqq 2 with M P N, the definition of AS and T S can be extended to any patch domain ω Ă Z 2 by replacing u by 9 u, its periodic extension to Z 2 . A similar extension can be derived in the finite continuous setting, i.e. Ω " T 2 .
Suppose we evaluate the scalar product auto-similarity AS sc pU, t, ωq with U a random field. Then the autosimilarity function is a random variable and its expectation depends on the second-order statistics of U . In the template case, the expectation of T S sc pU, v, t, ωq depends on the first-order statistics of U . This shows that auto-similarity and template similarity can exhibit very different behaviors even for the same similarity functions.
In the discrete case, it is well-known that, due to the curse of dimensionality, the L 2 norm does not behave well in large-dimensional spaces and is a poor measure of structure. Thus, considering u and v two images, s 2 pu, vq, the L 2 template similarity on full images, does not yield interesting information about the perceptual differences between u and v. The template similarity T S 2 pu, v, 0, ωq avoids this effect by considering patches which reduces the dimension of the data (if the cardinality of ω, denoted |ω|, is small) and also allows for fast computation of similarity mappings, see Figure 1 for a comparison of the different similarity functions on a natural image.
We extract patches from images as follow. For each position in the image we consider a square ω centered around this position. This operation is called patch lifting. In Figure 2, we investigate the behavior of patch lifting on different Gaussian random fields. Roughly speaking, patches are said to be similar if they are clustered in the patch space. Using Principal Component Analysis we illustrate that patches are more scattered in Gaussian white noise than in the Gaussian random field U " f˚W (with periodic convolution, i.e. f˚W pxq " ř yPΩ W pyq 9 f px´yq where 9 f is the periodic extension of f to Z 2 ), where W is a Gaussian white noise over Ω (a finite discrete grid) and f is the indicator function of a rectangle non reduced to a single pixel. . Gaussian models and spatial redundancy. In this experiment we illustrate the notion of spatial redundancy in two models. In (A), we present a 64ˆ64 Gaussian white noise. (B) Shows an indicator function f . In (C), we present a realization of the Gaussian random field defined by f˚W (with periodic convolution) where W is a Gaussian white noise over Ω (domain of size 64ˆ64). Note that f was chosen so that the two Gaussian random fields (A) and (C) have the same gray-level distribution for each pixel. To each pixel position in (A) and (C) we associate the surrounding patch, with patch domain ω (of size 3ˆ3q. Hence, for each image (A) and (C) we obtain 64ˆ64 " 5096 vectors each of size 3ˆ3 " 9. These 9-dimensional vectors are projected in a 3-dimensional space using Principal Component Analysis. In the subfigure (D), we display the 20 vectors closest to 0 in each case: Gaussian white noise model (in blue) and the Gaussian random field (C) (in red). The radius of the blue, respectively red, sphere represents the maximal L 2 norm of these 20 vectors in the Gaussian white noise model, respectively in model (C). Since the radius of the blue sphere is larger than the red one the points are more scattered in the patch space of (A) than in the patch space of (B). This implies that there is more spatial redundancy in (C) than in (A), which is expected.
We continue this investigation in Figure 3 in which we present the closest patches (of size 10ˆ10), for the L 2 norm, in two Gaussian random fields U " f˚W (where the convolution is periodic) for different functions f called spots, [25]. The more regular f is, the more similar the patches are. Limit cases are f " 0 (all patches are constant and thus all the patches are similar) and f " δ 0 , i.e. U " W .
We introduce the notion of autocorrelation. Let f P L 2 pZ 2 q. We denote by Γ f the autocorrelation of f , i.e. Γ f " f˚f where for any x P Z 2 ,f pxq " f p´xq and define the associated random field to a square-integrable function f as the stationary Gaussian random field U such that for any x P Ω E rU pxqs " 0 and Γpxq " Γ f pxq.
In Figure 4, we compare the patch spaces of natural images and the ones of their associated random fields. Since the associated Gaussian random fields lose all global structures, most of the spatial information is discarded. This situation can be observed in the patch space. In the natural images, patches containing the same highly spatial information (such as a white diagonal) are close for the L 2 norm. In Gaussian random field since this highly spatial information is lost, close patches for the L 2 norm are not necessarily perceptually close.

Asymptotic results
In this section we aim at giving explicit asymptotic expressions for the probability distribution functions of the auto-similarity and the template similarity in both discrete and continuous settings. Using general versions of the law of large numbers and central limit theorems we will derive Gaussian asymptotic approximations.
Additional assumptions are required in the case of template matching since we use an exemplar input image v to compute T S i pU, v, t, ωq. Let v P R Ω , where Ω is R 2 or Z 2 . We denote by pv k q kPN the sequence of the restriction a Gaussian spot f and a realization of the Gaussian random field U " f˚W , where the convolution is periodic and W is a Gaussian white noise. The associated random field is smooth and isotropic. The random field U " f˚W associated with a rectangular plateau f is no longer smooth nor isotropic. Images are displayed on the right of their respective spot. For each setting (Gaussian spot or rectangular spot) we present 12 patches of size 15ˆ15. In each case the topleft patch is the top-left patch in the presented realization of the random field, shown in green. Following from the top to the bottom and from the left to the right are the closest patches in the patch space for the L 2 norm. We discard patches which are spatially too close (if ω 1 and ω 2 are two patch domains we impose sup x,y }x´y} 8 ě 10).
of v to ω k , extended to Z 2 (or R 2 ) by zero-padding, i.e. v k pxq " 0 for x R ω k . We suppose that lim kÑ`8 |ω k | " 8, where |ω k | is the Lebesgue measure, respectively the cardinality, of ω k if Ω " R 2 , respectively Ω " Z 2 . Note that the following assumptions are well-defined for both continuous and discrete settings. The following assumption ensures the existence of spatial moments of any order for the function v.
Assumption 3.2 (A3.2). For any m, n P N, there exist β m P Rzt0u and γ m,n P R Ω such that (a) lim kÑ`8 |ω k | 1{2´| ω k |´1 ş ω k v 2m k pxqdx´β m¯" 0; (b) for any K Ă Ω compact, lim kÑ`8 sup xPK |ω k |´1 ş yPω k v 2m k pyqv 2n k px`yqdy´γ m,n pxq| " 0. . Natural images and Gaussian random fields. In this experiment we present the same image, f , which was used in Figure 1 and the associated Gaussian random field U " f˚W , where the convolution is periodic and W is a Gaussian white noise. As in Figure 3 we present under each image the top-left patch (of size 15ˆ15 and shown in green in the original images) and its 11 closest matches for the 2 similarity. We discard patches which are spatially too close (if ω 1 and ω 2 are two patch domains we impose sup x,y }x´y} 8 ě 10). Note that if a structure is clearly identified in the real image (black and white diagonals) and is retrieved in every patch, it is not as clear in the Gaussian random field.
Note that in the case where Ω is discrete, the uniform convergence on compact sets introduced in (b) is equivalent to the pointwise convergence. There exists γ P R Ω with for any K Ă Ω compact, lim kÑ`8 sup xPK |ω k |´1 ş yPω k v k pyq v k px`yqdy´γpxq| " 0.

Discrete case
In the discrete case, we consider a random field U over Z 2 and compute local similarity measurements. The asymptotic approximation is obtained when the patch size grows to infinity. In Theorems 3.4 and 3.6 we obtain Gaussian asymptotic probability distribution in the auto-similarity case and in the template similarity case. In Propositions 3.5 and 3.7 we give explicit mean and variance for the Gaussian approximations. We recall that N pµ, σ 2 q is the probability distribution of a Gaussian real random variable with mean µ and variance σ 2 .
Theorem 3.4 (Discrete case -asymptotic auto-similarity results). Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets defined for any k P N by, ω k " 0, m k ˆ 0, n k . Let f P R Z 2 , f ‰ 0 with finite support, W a Gaussian white noise over Z 2 and U " f˚W . For i " tp, pp, pq, sc, cosu with p P p0,`8q there exist µ i , σ 2 i , real valued functions on Z 2 , and pα i,k q kPN a positive sequence such that for any t P Z 2 z t0u we get The asymptotics derived in Theorem 3.4 can be extended to vectors of autosimilarities, i.e. selecting pt j q jPt1...N u a finite number of shifts the results of Theorem 3.4 hold for the sequence ppAS i pU, t j , ω k qq jPt1...N u q kPN . Note that in Theorem 3.4 if t varies with k such that for any k P N, pω k`tk q X ω k " H then similar results can be obtained with the usual law of large numbers and central limit theorem since true independence hold.
Proof. The proof is divided into three parts. First we show 1 and 2 for i " p, p and extends the result to i " p. Then we show 1 and 2 for i " sc. Finally, we show 1 and 2 for i " cos.
1. Let p P p0,`8q, t P Z 2 z t0u and define V p,t for any x P Z 2 by, V p,t pxq " |U pxq´U px`tq| p . We remark that for any k P N we have We first notice that U is R-independent with R ą 0, see Lemma B.5. Since for any x P Z 2 , V p,t pxq depends only on U pxq and U px`tq we have that V p,t is R t " R`}t} 8 -independent. Since U is stationary, so is V p,t . The random field V p,t admits moments of every order since it is the pth power of the absolute value of a Gaussian random field. Thus V p,t is a R t -independent second-order stationary random field. We can apply Lemma B.6 and we get (a) lim kÑ`8 |ω k | AS p,p pU, t, ω k q´µ p,p ptq¯" L N`0, σ 2 p,p ptq˘. with µ p,p ptq " E rV p,t p0qs and σ p,p ptq 2 " ř xPZ 2 Cov rV p,t pxq, V p,t p0qs. By continuity of the p-th root over r0,`8q we get 1 for i " p with By Lemma B.7 we get that E " pU p0q´U ptqq 2 ‰ " 2pΓ f p0q´Γ f ptqq ą 0 thus µ p,p ptq " E rV p,t p0qs ą 0. Since the pth root is continuously differentiable on p0,`8q we can apply the Delta method, see [14], and we get 2 for i " p with 2. We now prove the theorem for i " sc. Let t P Z 2 z t0u and define V sc,t for any x P Z 2 , V sc,t pxq " U pxqU px`tq. We remark that for any k P N we have Since for any x P Z 2 , V sc,t pxq depends only on U pxq and U px`tq, we have that V sc,t is R t " R`}t} 8independent. Since U is stationary, so is V sc,t . The random field V sc,t admits moments of every order since it is a product of Gaussian random fields. Thus V sc,t is a R t -independent second-order stationary random field. We can again apply Lemma B.6 and we get (a) lim kÑ`8 |ω k | AS sc pU, t, ω k q´µ sc ptq¯" L N`0, σ 2 sc ptq˘, with µ sc ptq " E rV sc,t p0qs and σ sc ptq 2 " ř xPZ 2 Cov rV sc,t pxq, V sc,t p0qs, which concludes the proof. 3. Finally, we consider the case i " cos. Let t P Z 2 z t0u and define V cos,t for any x P Z 2 , We remark that for any k P N we have The random field V cos,t admits moments of every order since it is a vector of products of Gaussian random fields. Thus V cos,t is a R t -independent second-order stationary random field. We can apply Lemma B.6 and there existμ cos ptq andC cos ptq such that (a) lim kÑ`8 a.s.μ cos ptq; We conclude the proof using the multivariate Delta method, [14].
In the following proposition we give explicit values for the constants involved in the law of large numbers and the central limit theorem derived in Theorem 3.4. We introduce the following quantities for k, P N and j P 0, k^ , where k^ " minpk, q, We also denote r j, " r j, , . Note that for all P N, r 0, " q 2 and ř j"0 r j, " q 2 . We also introduce the following functions: Note that ∆ f is a second-order statistic on the Gaussian field U " f˚W with W a Gaussian white noise over whereas r ∆ f is a fourth-order statistic on the same random field.
Proposition 3.5 (Explicit constants -Auto-similarity). In Theorem 3.4 we have the following constants for any t P Z 2 z t0u.
(i) If i " p with p " 2 and P N, then for all k P N, we get that α p,k " |ω k | 1{p2 q and where pr i,jk q i,j,kPN and pq k q kPN are given in (3.4). (ii) If i " sc, then for all k P N, we get that α sc,k " |ω k | and µ sc ptq " Γ f ptq and σ sc ptq 2 " (iii) if i " cos, then for all k P N, we get that α cos,k " 1 and Proof. The proof is postponed to Appendix D.
For example we have We now derive similar asymptotic properties in the template similarity case.
Theorem 3.6 (Discrete case -asymptotic template similarity results). Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets defined for any k P N, ω k " 0, m k ˆ 0, n k . Let f P R Z 2 , f ‰ 0 with finite support, W a Gaussian white noise over Z 2 , U " f˚W and let v, a real valued function on Z 2 . For i " tp, pp, pq, sc, cosu with p " 2 and P N, if i " p or pp, pq assume (A3.1) and (A3.2), if i " sc assume (A3.1) and (A3.3) and if i " cos assume (A3.1), (A3.2) and (A3.3). Then there exist µ i , σ 2 i P R and pα i,k q kPN a positive sequence such that for any t P Z 2 we get Note that contrarily to Theorem 3.4 we could not obtain such a result for all p P p0,`8q but only for even integers. Indeed, in the general case the convergence of the sequence`|ω k |´1E rT S p,p pU, v, t, ω k qs˘k PN , which is needed in order to apply Theorem B.3, is not trivial. Assuming that v is bounded it is easy to show that |ω k |´1E rT S p,p pU, v, t, ω k qs˘k PN is also bounded and we can deduce the existence of a convergent subsequence. In the general case, for Theorem 3.6 to hold with any p P p0,`8q, we must verify that for any t P Ω, there exist µ p,p ptq ą 0 and σ 2 p,p ptq ě 0 such that p,p ptq. We now turn to the proof of Theorem 3.6.
Proof. As for the proof of Theorem 3.4, the proof is divided into three parts. First we show 1 and 2 for i " pp, pq and extends the result to i " p. Then we show 1 and 2 for i " sc. Finally, we show 1 and 2 for i " cos.
1. Let p P p0,`8q, t P Z 2 and define V p,t the random field on Z 2 for any x P Z 2 , by V p,t pxq " |vpxq´U px`tq| p . We remark that for any k P N we have ‰ uniformly almost surely dominates V p,t pxq´E rV p,t pxqs. The random field V 8 p,t admits moments of every order since it is the pth power of the absolute value of a Gaussian random field and is stationary because U is. Thus V p,t is a R t -independent random field and V p,t pxq´E rV p,t pxqs is uniformly stochastically dominated by V 8 ‰ , a second-order stationary random field. Using (A3.2) and Lemma B.8, we can apply Theorem B.3 and B.4 and we get (a) lim kÑ`8 Note that since U is stationary we have for any t P Z 2 , µ p,p " µ p,p p0q " µ p,p ptq and σ 2 p,p " σ 2 p,p p0q " σ 2 p,p ptq. By continuity of the pth root over r0,`8q we get 1 for i " p with By Lemma B.8, we have that µ p,p ą 0. Since the pth root is continuously differentiable on p0,`8q we can apply the Delta method and we get 2 for i " p with 2. We now prove the theorem for i " sc. Let t P Z 2 and define V sc,t the random field on Z 2 such that for any x P Z 2 , V sc,t pxq "´vpxqU px`tq. We remark that for any k P N we have It is clear that for any k P N, T S sc pU, v, t, ω k q is a R-independent Gaussian random variable with E rT S sc pU, v, t, ω k qs " 0 and where we recall that v k is the restriction of v to ω k . The last sum is finite since Supp pf q finite implies that Supp pΓ f q is finite. Using (A3.3) we obtain that for any k P N, with µ sc " 0 and σ 2 sc " Since V sc,t is a R-independent second-order random field using (3.8) we can apply Theorems B.3 and B.4 to conclude. 3. We now consider the case i " cos. First, notice that In addition, using Slutsky's theorem and the fact that Proposition 3.7 (Explicit constants -template similarity). In Theorem 3.6 we have the following constants for any t P Z 2 .
(i) If i " p with p " 2 and P N, then we get that α p,k " |ω k | 1 p , and where pβ j q jPN , pγ i,j q i,jPN are given in (A3.2) and pr i,jk q i,j,kPN and pq k q kPN are given in (3.4). (ii) If i " sc then for all k P N, we get that α sc,k " |ω k | and µ sc " 0 , σ 2 sc " xγ, Γ f y.
(iii) If i " s cos then for all k P N, we get that α scos,k " 1 and Proof. The proof is postponed to Appendix D.
For example we have (3.10) Note that the limit mean and standard deviation do not depend on the offset anymore. Indeed, template similarity functions are stationary in t. If v has finite support then (A3.2) holds with β i " 0 and γ i,j " 0 as soon as i ‰ 0 or j ‰ 0. Remarking that β 0 " 1 and γ 0,0 " 1 we obtain that

Continuous case
We now turn to the continuous setting. Theorem 3.8, respectively Theorem 3.10, is the continuous counterpart of Theorem 3.4, respectively Theorem 3.6.
Theorem 3.8 (Continuous case -asymptotic auto-similarity results). Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets defined for any k P N by, ω k " r0, m k sr 0, n k s. Let U be a zero-mean Gaussian random field over R 2 with covariance function Γ. Assume (A2.2) and that Γ has finite support. For i P tp, pp, pq, sc, cosu with p P p0,`8q there exist µ i , σ 2 i , real valued functions on R 2 , and pα i,k q kPN a positive sequence such that for any t P R 2 z t0u we get Proof. The proof is the same as the one of Theorem 3.4 replacing Lemmas B.6 and B.7 by Lemma C.3 and Lemma C.4.
Proposition 3.9 (Explicit constants -Continuous auto-similarity). Constants given in Proposition 3.5 apply to Theorem 3.8 provided that Γ f is replaced by Γ in (3.5).
Proof. The proof is the same as the one of Proposition 3.5.
Theorem 3.10 (Continuous case -asymptotic template similarity results). Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets defined for any k P N by, ω k " r0, m k sr 0, n k s. Let U be a zero-mean Gaussian random field over R 2 with covariance function Γ. Assume (A2.2) and that Γ has finite support. For i P tp, pp, pq, sc, cosu with p P p0,`8q, if i " p or pp, pq assume (A3.1) and (A3.2), if i " sc assume (A3.1) and (A3.3) and if i " cos assume (A3.1), (A3.2) and (A3.3). Then there exist µ i , σ 2 i P R and pα i,k q kPN a positive sequence such that for any t P R 2 we get Proof. The proof is the same as the one of Theorem 3.6. Proposition 3.11 (Explicit constants -Continuous auto-similarity). Constants given in Proposition 3.7 apply to Theorem 3.10 provided that Γ f is replaced by Γ in (3.5).
Proof. The proof is similar to the one of Proposition 3.7.

Speed of convergence
In the discrete setting, Theorem 3.4 justifies the use of a Gaussian approximation to compute AS i pU, t, ωq. However this asymptotic behavior strongly relies on the increasing size of the patch domains. We define the patch size to be |ω|, the cardinality of ω, and the spot size |Supp pf q | to be the cardinality of the support of the spot f . The quantity of interest is the ratio r " patch size spot size . If r " 1 then the Gaussian random field associated to f can be well approximated by a Gaussian white noise from the patch perspective. If r « 1 this approximation is not valid and the Gaussian approximation is no longer accurate, see Figure 5. We say that an offset t is detected in a Gaussian random field if AS i pU, t, ωq ď aptq for some threshold aptq. In the experiments presented in Figure 6  Scalar product auto-similarities and squared L 2 auto-similarities are computed for a fixed offset p70, 100q. We then plot the normalized histogram of these values. The red curve corresponds to the standard Gaussian N p0, 1q. On the top row r " 100 " 1 and the Gaussian approximation is valid. On the bottom row r « 1 and the Gaussian approximation is not valid.
and Table 1 the threshold is given by the asymptotic Gaussian inverse cumulative distribution function evaluated at some quantile. The parameters of the Gaussian random variable are given by Proposition 3.5. We find that except for small spot sizes and large patches, i.e. r " 1, the approximation is not valid. More precisely, let U " f˚W with f a finitely supported function over Z 2 and W a Gaussian white noise over Z 2 . Let ω Ă Z 2 and let Ω 0 be a finite subset of Z 2 . We compute ř tPΩ0 1 ASipU,t,ωqďaptq , with aptq defined by the inverse cumulative distribution function of quantile 10{|Ω 0 | for the Gaussian N pµ, σ 2 q where µ, σ 2 are given by Theorem 3.4 and Proposition 3.5. Note that aptq would satisfy P rAS i pU, t, ωq ď aptqs « 10{|Ω 0 | if the approximation for the cumulative distribution function was correct. In other words, if the Gaussian asymptotic was always valid, we would have a number of detections equal to 10 independently of r. This is clearly not the case in Table 1. One way to interpret this is by looking at the left tail of the approximated distribution for s 2,2 and s sc on Figure 5. For s sc the histogram is above the estimated curve, see (a) in Figure 6 for example. Whereas for s 2,2 the histogram is under the estimated curve. Thus for s sc we expect to obtain more detections than what is predicted whereas we will observe the opposite behavior for s 2,2 . This situation is also illustrated for similarities Figure 6. Theoretical and empirical cumulative distribution Function. This experiment illustrates the non-Gaussianity in Figure 5. In both cases, the red curve is the inverse cumulative distribution function of the standard Gaussian and the blue curve is the empirical inverse cumulative distribution function of normalized auto-similarity functions computed with 10 4 samples of Gaussian models. We present auto-similarity results obtained for t " p70, 100q and similarity function s sc (on the left) and s 2 (on the right). We note that for rare events, see the magnified region, the theoretical inverse cumulative distribution function is above the empirical inverse cumulative distribution function. The opposite behavior is observed for similarity s 2 . These observations are in accordance with the findings of Table 1.  ) auto-similarity function. We only consider patch domains larger than spot domains. We generate 5000 Gaussian random field images of size 256ˆ256 for each setting (with spot the indicator of the spot domain). We set α " 10{256 2 . For each setting we compute aptq the inverse cumulative distribution function of N pµ i ptq, σ 2 i ptqq evaluated at quantile α, with µ i and σ 2 i given by Proposition 3.5. For each pair of patch size and spot size we compute ř tPΩ 1 ASipu,t,ωqďaptq , namely the number of detections, for all the 5000 random fields samples. The empirical averages are displayed in the table. If AS i pu, t, ωq had Gaussian distribution with parameters given by Proposition 3.5 then the number in each cell would be ř tPΩ P rAS i pU, t, ωq ď aptqs « 10. s 2 and s sc in Figure 6 in which we compare the asymptotic cumulative distribution function with the empirical one.
In the next section we address this problem by studying non-asymptotic cases for the s 2,2 auto-similarity function in both continuous and discrete settings.

Discrete periodic case
In this section Ω is a finite rectangular domain in Z 2 . We fix ω Ă Ω. We also define f a function over Ω. We consider the Gaussian random field U " f˚W (we consider the periodic convolution) with W a Gaussian white noise over Ω.
In the previous section, we derived asymptotic properties for similarity functions. However, a necessary condition for the asymptotic Gaussian approximation to be valid is for the spot size to be very small when compared to the patch size. This condition is not often met and non-asymptotic techniques must be developed. For instance it should be noted that the distribution of the s sc template similarity, T S sc pU, v, t, ωq, is Gaussian for every ω. We might also derive a non-asymptotic expression for the template similarity in the cosine case if the Gaussian model is a white noise model. In what follows we restrict ourselves to the auto-similarity framework and consider the square of the L 2 norm auto-similarity function, i.e. AS 2,2 pU, t, ωq. In this case we present an efficient method to compute the cumulative distribution function of the auto-similarity function even in the non-asymptotic case.
Proposition 4.1 (Squared L 2 auto-similarity function exact probability distribution function). Let Ω " pZ{M Zq 2 with M P N, ω Ă Ω, f P R Ω and U " f˚W where W is a Gaussian white noise over Ω. The following equality holds for any t P Ω up to a change of the underlying probability space with Z k independent chi-square random variables with parameter 1 and λ k pt, ωq the eigenvalues of the covariance matrix C t associated with function ∆ f pt,¨q, see equation (3.5), restricted to ω, i.e. for any x 1 , x 2 P ω, C t px 1 , x 2 q " ∆ f pt, x 1´x2 q.
Proof. Let t P Ω and V t be given for any x P Ω by V t pxq " U pxq´U px`tq. It is a Gaussian vector with mean 0 and covariance matrix C V given for any x 1 , x 2 P Ω by The covariance of the random field P ω pV t q, the restriction of V t to ω, is given by the restriction of C V to ω. This new covariance matrix, C t , is symmetric and the spectral theorem ensures that there exists an orthonormal basis B such that C t is diagonal when expressed in B. Thus we obtain that P ω pV t q " ř e k PB xP ω pV t q, e k ye k . It is clear that, for any k P 0, |ω|´1 , xP ω pV t q, e k y is a Gaussian random variable with mean 0 and variance e T k C t e k " λ k pt, ωq ě 0. We set K " tk P 0, |ω|´1 , λ k pt, ωq ‰ 0u and define X a random vector in R |ω| such that where X K´i s the restriction of X to the indices of K´" 0, |ω|´1 zK and Y is a standard Gaussian random vector on R |K´| independent from the sigma field generated by tpX k q, k P Ku. By construction we have E rX k X s " 0 if P K and k P K´, or P K´and k P K´. Suppose now that k, P K. We obtain that

. Blue circles correspond to the 100 eigenvalues computed with
Matlab routine for offset p5, 5q in (B), respectively p10, 20q in (C), and red crosses correspond to the 100 approximated eigenvalues for the same offsets. Note that a standard routine takes 273s for 10ˆ10 patches on 256ˆ256 images whereas it only takes 1.11s when approximating the eigenvalues using the discrete Fourier transform.
Setting Z k " X 2 k concludes the proof.
Note that if ω " Ω then we obtain that the covariance matrix C t is block-circulant with circulant blocks and the eigenvalues are given by the discrete Fourier transform.
In order to compute the true cumulative distribution function of the auto-similarity square L 2 norm we need to: (1) compute the eigenvalues of a covariance matrix in M |ω| pRq; (2) compute the cumulative distribution function of a positive-weighted sum of independent chi-square random variable with weights given by the computed eigenvalues. Storing all covariance matrices for each offset t is not feasible. For instance considering a patch of size 10ˆ10 and an image of size 512ˆ512 we have approximately 2.6ˆ10 9 coefficients to store, i.e. 10.5GB in float precision. In the rest of the section we suppose that t and ω are fixed and we denote by C t the covariance matrix associated to the restriction of ∆ f pt,¨q to ω`p´ωq. In Proposition 4.2 we propose a method to approximate the eigenvalues of C t by using its specific structure. Indeed, as a covariance matrix, C t is symmetric and positive and, since its associated Gaussian random field is stationary, it is block-Toeplitz with Toeplitz blocks, i.e. is block-diagonally constant and each block has constant diagonals. In the one-dimensional case these properties translate into symmetry, positivity and Toeplitz properties of the covariance matrix. Proposition 4.2 is stated in the one-dimensional case for the sake of simplicity but two-dimensional analogous can be derived. Note that this approximation is not always sharp as shown in Figure 7.
We recall that the Frobenius norm of a matrix of size nˆn is the L 2 norm of the associated vector of size n 2 .
Proposition 4.2 (Eigenvalues approximation). Let b be a function defined over ´pn´1q, n´1 with n P Nz t0u. We define T b pj, q " bpj´ q for j, P 0, n´1 . The matrix T b is a circulant matrix if and only if b is n-periodic. T b is symmetric if and only if b is symmetric. Let b be symmetric, defining ΠpT b q the projection of T b onto the set of symmetric circulant matrix for the Frobenius product, we obtain that 1. the projection satisfies ΠpT b q " T c with cpjq "`1´j n˘b pjq`j n bpn´jq for all j P 0, n´1 and c is extended by n-periodicity to Z; 2. the eigenvalues of ΠpT b q are given by´2 Repdpjqq´bp0q¯j P 0,n´1 with dpjq "`1´j n˘b pjq, andd is the discrete Fourier transform over 0, n´1 ; 3. let pλ j q jP 1,n be the sorted eigenvalues of T b and pλ j q jP 1,n the sorted eigenvalues of ΠpT b q (in the same order). For any j P 1, n , we have |λ j´λj | ď }T b´Π pT b q} Fr ; 4. if T b is positive-definitive then ΠpT b q is positive-definite.
Proof. (1) Let T c be an element of the symmetric circulant matrices set. Minimizing }T b´Tc } 2 Fr in cpjq jP 0,n´1 we get that cpjq satisfies for any j P 0, n´1 cpjq " argmin sPR`2 pn´jqps´bpjqq 2`2 jps´bpn´jqq 2˘, which gives the result.
(2) Since T c " ΠpT b q is circulant, its eigenvalues are given by the discrete Fourier transform of c. We have that if i ‰ 0 then cpiq " 9 dpjq`9 dp´jq with dpjq "`1´j n˘b pjq and 9 d its extension to Z by n-periodicity. We also have cp0q " bp0q. We conclude the proof by taking the discrete Fourier transform of c.
(3) The proof of the Lipschitz property on the sorted eigenvalues of symmetric matrices with respect to the L 2 matricial norm can be found in [13]. We conclude using the fact that the L 2 matricial norm is upper-bounded by the Frobenius norm.
(4) This result is a special case of the spectrum contraction property of the projection proved in Theorem 2 of [11].
In Figure 7 we display the behavior of the projection for the eigenvalues in the two-dimensional case. The measure we consider is the Normalized Root Mean Square Deviation withλ k pt, ωq the approximation of the eigenvalues, for every possible offset in the image and λ k pt, ωq the true eigenvalues, for every possible offset. Computing the eigenvalues of the projection is done via Fast Fourier Transform (FFT) which is faster than standard routines for computing eigenvalues of Toeplitz matrices. The major cons of using such an approximation is that it may not be valid for small offsets t P Ω as shown in Figure 7. However, in most cases the random field is smooth and in this case, see Figure 8, the approximation is satisfactory. We also highlight that for similarity detection purposes, see Figure 9, the level of precision achieved by our approximation is satisfactory, see [7]. Suppose the approximation of the eigenvalues is valid, we need an efficient algorithm to compute the distribution of the associated positive-weighted sum of chi-square random variables in equation (4.1). Exact computation has been derived by Imhof in [31] but requires to compute heavy integrals. This exact method, named Imhof method in the following, will be used as a baseline for other algorithms. Numerous methods such as differential Figure 8. Eigenvalues approximation. Same study as the one conducted in Figure 7 with f " 1 t1,2,3u 2 . Note that in this case the approximation is better than the one presented in Figure 7. Figure 9. Similarity detection. In this figure we illustrate the accuracy of the different proposed approximations of the cumulative distribution function of AS 2,2 pU, t, ωq. We say that an offset t is detected in an image if AS 2,2 pu, t, ωq ď aptq for some threshold aptq P R. In every image, in green we display the patch domain ω (in the center of the image) and in red we display the shifted patch domain for detected offsets with function aptq such that for any t P Ω, P rAS 2,2 pU, t, ωq ď aptqs " 1{256 2 , where U is given by the Gaussian random field f˚W where f is the original image of fabric and W is a Gaussian white noise over Ω " 256ˆ256. Approximations of the cumulative distribution function of AS 2,2 pU, t, ωq lead to approximations of aptq. The most precise approximation is given in (A) where the eigenvalues are computed using a Matlab routine and the cumulative distribution function is given by the Imhof method. In (B) we approximate the eigenvalues using the projection described in Proposition 4.2 and still use the Imhof method. It yields twice as many detections. In (C) Wood F method is used instead of Imhof's yielding less detections but performing seven times faster. Interestingly errors seem to compensate and the obtained result with Wood F method is very close to the results obtained with the baseline algorithm in (A). In (D) HBE method is used instead of Imhof's, in this case we obtain too many detections, i.e. the approximation of the cumulative distribution function is not valid.
equations [19], series truncation [38], negative binomial mixtures [46] approaches were later introduced but all require stopping criteria such as truncation criteria which can be hard to set. We focus on cumulant methods which generalize and refine the Gaussian approximations used in Section 3. These methods rely on computing moments of the original distribution and then fitting a known probability distribution function to the objective distribution using these moments. Bodenham et al. in [6] show that the following methods can be efficiently computed: -Gaussian approximation (discarded due to its poor results for small patches as illustrated in Sect. 3); -Hall-Buckley-Eagleson [10,29] (HBE), (three moments fitted Gamma distribution); -Wood F [56] (three moments fitted Fischer-Snedecor distribution).
Other methods such as the Lindsay-Pilla-Basak-4 method, which relies on the computation of eight moments, are slower than HBE by a factor 350 at least, see [6], and are thus discarded. In Figure 9 we investigate the trade-off between computational speed and accuracy of these methods for the task of detection.
The experiments conducted in Figure 9 show that the HBE approximation does not give good results when evaluating the probability of rare events. This was already noticed by Bodenham et al. in [6] who stated that "Hall-Buckley-Eagleson method is recommended for most practitioners [. . .]. However, [. . .], for very small probability values, either the Wood F or the Lindsay-Pilla-Basak method should be used".

Continuous periodic case
To conclude we show that a similar non-asymptotic study can be conducted in continuous settings. Let Ω " T 2 , ω Ă Ω and let U be a zero-mean Gaussian random field on Ω with covariance function Γ. Assume (A2.2), then the following equality holds for any t P Ω up to a change of the underlying probability space with Z k independent chi-square random variables with parameter 1 and λ k pt, ωq the eigenvalues of the kernel C t associated with function ∆pt,¨q " 2Γptq´Γp¨`tq´Γp¨´tq restricted to ω, i.e. for any x 1 , x 2 P ω, C t px 1 , x 2 q " ∆pt, x 1´x2 q.
Proof. We consider the stationary Gaussian random field P ω pV t q over ω defined by the restriction to ω where for any x P Ω by V t pxq " U pxq´U px`tq. The Karhunen-Loeve theorem [28] ensures the existence of pλ k pt, ωqq kPN P R Ǹ , pX k q kPN a sequence of independent normal Gaussian random variables and pe k q kPN a sequence of orthonormal function over L 2 pωq such that We define the sequence pI n q nPN " p ş ω´ř n k"0 a λ k pt, ωqe k pxqX k¯2 dxq nPN . We have, using the Cauchy-Schwarz inequality on L 2 pAˆωq and (4.3) where we used the Fubini theorem in the last inequality. Using the dominated convergence theorem in (4.4) with integral domination given by sup nPN sup xPω E " pP ω pV t qpxq´ř n k"0 a λ k pt, ωqe k pxqX k q 2 ı we conclude that pI n q nPN converges to AS 2,2 pU, t, ωq in L 1 pAq. Thus there exists a subsequence of pI n q nPN which converges almost surely to AS 2,2 pU, t, ωq. We also have I n " ş ω´ř n k"0 a λ k pt, ωqe k pxqX k¯2 dx " ř n k"0 λ k pω, kqX 2 k by orthonormality and thus the sequence pI n q nPN is almost surely non-decreasing. We get that pI n q nPN converges almost surely to AS 2,2 pU, t, ωq which can be rewritten as The characterization of pλ k pt, ωq, e k pxqq is given by the Karhunen-Loeve theorem and e k pxq is solution of the following Fredholm equation for all x P ω ż ω ∆pt, x´yqe k pyq dy " λ k pt, ωqe k pxq.
Setting Z k " X 2 k concludes the proof. Note that if ω " T 2 then the solution of the Fredholm equation is given by the Fourier series of Γ.

Appendix A. Multidimensional central limit theorems
In this section we provide an extension of ( [35], Thm. 2) to the multidimensional case. We recall the notion of dependency graph as introduced in [35]. Let pX i q iPN be R d -valued random variables. A graph is a dependency graph for pX i q iPN if the two following conditions are satisfied.
1. There is a one-to-one correspondence between pX i q iPN and the vertices of the graph. 2. If two sets of vertices are not connected then the corresponding random variables are independent.
Theorem A.1. Let pX i,j q pi,jqPN 2 be a sequence of R d -valued random variables and pN n q nPN P N N . For any n P N, assume that there exists A n , M n ě 0 such that for any j P N, }X n,j } ď A n and that the dependency graph of pX n,j q jPN is of degree M n at most. For any n P N let S n " ř Nn j"1 X n,j and C n " Cov rS n s. Assume that there exists m 0 P N and C P M d pRq such that for any n P N, lim nÑ`8 pN n {M n q 1{m0 M n A n " 0 and lim nÑ`8 C n " C. Then, S n´E rS n s converges (in the weak sense) towards N p0, Cq.
Proof. Let a P R d and consider pX a i,j q pi,jqPN 2 such that for any i, j P N, X a i,j " xX i,j , ay. We also introduce for any n P N, S a n " ř Nn j"1 x a n,j . Assume that a J Ca " 0. Then, using the Bienaymé -Tchebychev inequality, we have for any ε ą 0 lim nÑ`8 P r|S a n´E rS a n s | ą εs ď lim Hence, xa, S n´E rS n sy converges (in the weak sense) towards xa, Zy with Z a d-dimensional Gaussian random variable with zero mean and covariance matrix C. If a J C n a ‰ 0 then using [35] we have that xa, S n´E rS n sy converges (in the weak sense) towards xa, Zy. We conclude using the Cramér-Wold theorem, [15], Theorem 1.
Similarly to ( [35], Thm. 2), we can replace the condition }X n,i } ď A n by the following condition: for any a P R d with a ‰ 0 Indeed, this implies that for any a P R d , lim nÑ`8 M n ř Nn j"1 E " xa, X n,j y 2 1 |xa,Xn,j y|ąAn ‰ " 0, which is the Lindeberg type condition identified in ( [35], Thm. 2).

Appendix B. Asymptotic theorems -discrete case
We start by introducing two notions which will be crucial in order to derive a law of large numbers and a central limit theorem in broad settings. The R-independence, see Definition B.1, ensures long-range independence whereas stochastic domination will replace integrability conditions in the standard law of large numbers or central limit theorem.
The notion of R-independence generalizes to R 2 and Z 2 the associated one-dimensional concept, see [5] and its extension to N 2 [44,53].
Definition B.1 (R-independence). Let d P N, Ω " R 2 or Ω " Z 2 and V be a d-dimensional random field over Ω. Let K 1 , K 2 Ă Ω be two compact sets, and V | Ki be the restriction of V to K i , i P t1, 2u. We say that V is R-independent, with R ě 0, if V | K1 is independent from V | K2 as soon as d 8 pK 1 , K 2 q " min xPK1,yPK2 Note that in the case of Ω " Z 2 , compacts sets K 1 and K 2 are finite sets of indices. This notion of Rindependence will replace the traditional assumption of independence in asymptotic theorems.
Definition B.2 (Uniform domination). Let Ω " R 2 or Ω " Z 2 and let V, r V be real random fields over Ω. We say that: (a) r V uniformly stochastically dominates V if for any α ě 0 and x P Ω, P rV pxq ě αs ď P " r V pxq ě α ı ; (b) r V uniformly almost surely dominates V if for any x P Ω, V pxq ď r V pxq almost surely.
Note that if r V uniformly almost surely dominates V then r V uniformly stochastically dominates r V . The following theorem is a two-dimensional law of large numbers with weak dependence assumptions. It is a slight modification of Corollary 4.1 (ii) in [53].
Theorem B.3. Let d P N. Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets such that for any k P N, ω k " 0, m k ˆ 0, n k . Let V be a d-dimensional R-independent random field over Z 2 , with R ě 0, such that }V pxq´E rV pxqs } is uniformly stochastically dominated by r V , a real second-order stationary random field over Z 2 . Then V is a second-order random field. In addition, assume that there exists µ P R d such that lim kÑ`8 |ω k |´1 ř xPω k E rV pxqs " µ. Then it holds that Proof. Without loss of generality we can suppose that d " 1 and that for any x P Z 2 , E rV pxqs " 0. In order to apply Corollary 4.1 (ii) in [53] we must check that: (a) V is R´independent; (b) |V | is uniformly stochastically dominated by a random field r V and there exists r P r1, 2r such that for any Item (a) is given in the statement of Theorem B.3 and |V | is uniformly stochastically dominated by the random field r V 0 defined for any x P Z 2 by r Using that lim kÑ`8 |ω k |´1 ř xPω k E rU pxqs " µ, we conclude the proof.
We now turn to an extension of the central limit theorem to two-dimensional random fields with weak dependence assumptions. This result is a consequence of Theorem A.1.
Theorem B.4. Under the hypotheses of Theorem B.3 and assuming that there exist µ P R d and C P M d pRq such that Then it holds that Proof. For any i, j P N, let X i,j " pV px j q´E rV px j qsq|ω i |´1 {2 with px j q jPN such that for any k P N, tV px j q, j P 1, |ω k | u " tV pxq, x P ω k u. For any n P N, let N n " |ω n |. Then, we have that for any n P N, ř Nn j"1 X n,j " |ω k |´1 {2 ř xPω k pV pxq´E rV pxqsq. Since V is R-independent each vertex of the dependency graph of pX i,j q i,jPN 2 has its degree bounded by p2R`1q 2 and therefore for any n P N, M n " p2R`1q 2 . For any n P N, let A n " |ω n | α with α P p1{3, 1{2q. Using thatṼ uniformly stochastically dominates p}V pxq´E rV pxqs }q xPZ 2 we obtain that for any a P R d ı .
Hence, since lim nÑ`8 A 2 n |ω n | " lim nÑ`8 |ω n | 1´2α "`8 we get that lim nÑ`8 Letting m 0 " 3 we get that In addition, we have that for any n P N C n " Cov The following lemma explicits a class of Gaussian random fields over Z 2 such that the R-independence property holds for some R ě 0.
Lemma B.5. Let f P R Z 2 with finite support Supp pf q Ă ´r, r 2 , where r P N. Let W be a Gaussian white noise over Z 2 and V " f˚W then V is a R-independent second-order random field with R " 2r.
Proof. V is a Gaussian random field such that for any x, y P Z 2 E rV pxqs " 0 , Cov rV pxq, V pyqs " Note that since Supp pf q Ă ´r, r we have Supp pΓ f q Ă ´R, R with R " 2r. For any x, y P Z 2 such that }x´y} 8 ą R, using (B.7), we obtain Cov rV pxq, V pyqs " Γ f px´yq " 0. (B.8) Let K 1 , K 2 Ă Z 2 two finite sets with sup xPK1,yPK2 }x´y} 8 ą R and consider V | Ki the restriction of V to K i for i " t1, 2u. Using (B.8), we get that for any x P K 1 , y P K 2 we have Cov rV | K1 pxq, V | K2 pyqs " 0.
The following lemma gives specific conditions on random fields in order for Theorems B.3 and B.4 to hold.
Lemma B.6. Let d P N. Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets given for any k P N by, ω k " 0, m k ˆ 0, n k . Let V be a d-dimensional R-independent second-order stationary random field over Z 2 , with R ě 0. Then for all k P N (a) |ω k |´1 ř xPω k E rV pxqs " E rV p0qs ; (b) lim kÑ`8 |ω k |´1 ř x,yPω k Cov rV pxq, V pyqs " ř xPZ 2 Cov rV pxq, V p0qs . In addition, equations (B.1) and (B.2) hold with µ " E rV p0qs and C " ř xPZ 2 Cov rV pxq, V p0qs which is finite. Proof. Item (a) is immediate by stationarity. Concerning (b), for any k P N we have by stationarity where g k P R Z 2 satisfies for any x P Z 2 , g k pxq " |ω k |´11 ω k˚1 ω k pxq. For any k P N, x P Z 2 we have 0 ď g k pxq ď 1 and lim kÑ`8 g k pxq " 1. For any x P Z 2 such that }x} 8 ą R, Cov rV pxq, V p0qs " 0 and then ř xPZ 2 | Cov rV pxq, V p0qs | ă`8. Proof. For any t P Z 2 , let τ t f " f p¨`tq. By the definition of the autocorrelation Γ f and using the Cauchy-Schwarz inequality we get that for any t P Z 2 with equality if and only if f " ατ t f , with α ‰ 0 since f ‰ 0. This implies that Supp pτ t pf qq " Supp pf q. As a consequence t " 0, which concludes the proof.
The following lemma ensures that items (a) and (b) in Theorem B.4 are satisfied in the template similarity case when imposing summability conditions over v.
Note that this lemma is also valid in the continuous case.

Appendix C. Asymptotic theorems -continuous case
We now turn to the continuous setting. We start by stating the continuous counterparts of Theorems B.3 and B.4. The following theorem, given here for completeness, can be found with different assumptions (in the one-dimensional case) in [42].
Theorem C.1. Let d P N. Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets given for any k P N by, ω k " r0, m k sˆr0, n k s. Let V be a d-dimensional R-independent random field over R 2 , with R ě 0, such that }V pxq´E rV pxqs } is uniformly stochastically dominated by r V , a stationary random field of order r ą 2 over R 2 . Then V is a second-order random field. In addition, assume V is sample path continuous and that there exists µ P R d given by lim kÑ`8 |ω k |´1 ş xPω k E rV pxqs dx " µ. Then it holds that µ. (C.1) Proof. Without loss of generality we can suppose that d " 1 and that for any x P Ω, E rV pxqs " 0. Let pσ k q kPN P R N given for any k P N by with Ω k " r0, ks 2 . Since V is R-independent, for any x, y P Ω such that }x´y} 8 ą R, we have Cpx, yq " 0. Hence for k large enough we obtain Using that r V uniformly stochastically dominates |V |, the stationarity of r V , and the Cauchy-Schwarz inequality, we obtain for any x, y P Ω, sup a.s.
Theorem C.2. Under the hypotheses of Theorem C.1 and assuming that there exist µ P R d and C P M d pRq such that (a) lim kÑ`8 |ω k |´1 {2 ş xPω k pE rV s pxq´µq dx " 0 ; (b) lim kÑ`8 |ω k |´1 ş x,yPω k Cov rV pxq, V pyqs dxdy " C. Then it holds that Proof. Let a P R d . We consider the d-dimensional random field ξ over R 2 defined for any x P R 2 by ξpxq " V pxq´E rV pxqs. We define also the weight functions pg n q nPN given for any n P N by g n pxq " |ω n |´1 {2 1 xPωn . For any n P N, let S n " ş R 2 g n pxqξpxqdx. We have for any n P N, S n " |ω n |´1 {2 ż xPωn pV pxq´E rV pxqsqdx.
Let ξ a be the one-dimensional random field over R 2 such that for any xÃ P R 2 , ξ a pxq " xa, ξpxqy and pS a n q nPN be the sequence of real-valued random variables such that for any n P N, S a n " xa, S n y. Then for any n P N, S a n " ş R 2 g n pxqξ a pxqdx. Using (b) we have that lim nÑ`8 E " pS a n q 2 ‰ " lim By assumption, ξ a is stochastically dominated by }a} r V and therefore for any x P R 2 we have E r|ξ a | r s ă`8. We conclude the proof upon using the Cramér-Wold theorem ( [15], Thm. 1).
The following lemmas are the continuous versions of Lemma B.6 and B.7.
Lemma C.3. Let d P N. Let pm k q kPN , pn k q kPN be two positive increasing integer sequences and pω k q kPN be the sequence of subsets such that for any k P N, ω k " r0, m k sˆr0, n k s. Let V be a d-dimensional, R-independent random field of order r ą 2 over R 2 , with R ě 0. Assume that V is sample path continuous, then for all k P N (a) |ω k |´1 ş xPω k E rV pxqs dx " E rV p0qs ; (b) lim kÑ`8 |ω k |´1 ş x,yPω k Cov rV pxq, V pyqs dxdy " ş xPR 2 Cov rV pxq, V p0qs dx. In addition, (C.1) and (C.5) hold with µ " E rV p0qs and C " ş xPR 2 Cov rV pxq, V p0qs dx. Proof. (a) The proof is immediate since for any x P R 2 , E rV pxqs " E rV p0qs.
By the Fubini-Lebesgue theorem we obtain that for any k P N, Cov rV pxq, V pyqs dxdy " ż xPR 2 Cov rV pxq, V p0qs g k pxqdx , where g k P L 8 pR 2 q satisfies for any x P R 2 , g k pxq " |ω k |´11 ω k˚1 ω k pxq. For any k P N, x P R 2 we have 0 ď g k pxq ď 1 and lim kÑ`8 g k pxq " 1. For any x P R 2 such that }x} 8 ą R t , Cov rV pxq, V p0qs " 0 and then ż xPR 2 | Cov rV pxq, V p0qs |dx ă`8.
Using the dominated convergence theorem we get that |ω k |´1 ż x,yPω k Cov rV pxq, V pyqs dxdy " ż xPR 2 Cov rV pxq, V p0qs dx , Since V is R-independent we conclude the proof by applying Theorem C.1 and C.2.
Lemma C.4. Let Γ be a function over R 2 , Γ ‰ 0, such that for any x, y P R 2 , Cpx, yq " Γpx´yq with C the covariance function of V a second-order random field over R 2 . Assume that Γ has finite support. Then it holds for any t P R 2 , Γptq ď Γp0q, with equality if and only if t " 0.
Proof. Upon replacing for any x P R 2 , V pxq by V pxq´E rV pxqs we suppose that E rV pxqs " 0. Using the Cauchy-Schwarz inequality and the stationarity of V we get for any t P R 2 and x P R 2 Γptq " E rV px`tqV pxqs ď E " with equality if and only if V px`tq " αpxqV pxq with αpxq P R. Since V is stationary and V ‰ 0 we get that for any x, t P R 2 , E " V px`tq 2 ‰ " E " V pxq 2 ‰ ą 0. Thus αpxq 2 " 1 and for all n P N, V pntq "˘V p0q. If t ‰ 0 then there exists n P N such that nt R Supp pΓq and then we have 0 " Γpntq " E rV pntqV p0qs "˘E " V p0 2 q ‰ ‰ 0 , which is absurd. Thus the equality in the inequality holds if and only if t " 0.

Appendix D. Explicit constants
In order to derive precise constants in Theorems 3.4 and 3.6 we use the following lemma which is a consequence of the Isserlis formula [32].
2. Let i " sc, t P Z 2 z t0u and V sc,t be a Gaussian random field given for any x P Z 2 , by V t pxq " U pxqU px`tq.