Loïc Denis

André Ferrari

David Mary

Laurent Mugnier

Eric Thiébaut

Fast and robust detection of a known pattern in an image

Keywords: robust detection, locally most powerful test (LMP), Cauchy distribution

Many image processing applications require to detect a known pattern buried under noise. While maximum correlation can be implemented efficiently using fast Fourier transforms, detection criteria that are robust to the presence of outliers are typically slower by several orders of magnitude. We derive the general expression of a robust detection criterion based on the theory of locally optimal detectors. The expression of the criterion is attractive because it offers a fast implementation based on correlations. Application of this criterion to Cauchy likelihood gives good detection performance in the presence of outliers, as shown in our numerical experiments. Special attention is given to proper normalization of the criterion in order to account for truncation at the image borders and noise with a non-stationary dispersion.

I. INTRODUCTION

We consider the problem of detecting a known pattern from a noisy image. In its simplest formulation, this detection problem amounts to a binary hypothesis choice:

H 0 : y = H 1 : y = α m(x) + , with α > 0 (1)
where under H 0 the observation y ∈ R n corresponds to the realization of a random process accounting for noise. Under hypothesis H 1 , the observation is the superimposition of noise and the pattern m(x) ∈ R n of (unknown) amplitude α > 0 and location x.

In many practical cases, there are several observations y i that display strong deviations from the model α m i (x). These deviations may be due to non-linearities of the sensor (e.g., saturation, dead/hot pixel) or mis-modeling (e.g., structured background). The random vector should not then be modeled by a white Gaussian distribution, but rather by a nonstationary, non-Gaussian, and perhaps correlated distribution. In this work, we keep the simplifying independence assumption but consider non-Gaussian noise distributions with heavy tails in order to derive a robust detection method.

When the measurement and modeling errors are centered and can be considered Gaussian, with a known covariance matrix Γ, the generalized likelihood ratio test takes the form of a matched filter [START_REF] Kay | Fundamentals of statistical signal processing[END_REF]: [START_REF] Fitch | Fast robust correlation[END_REF] when taking into account the prior information that α ≥ 0. Evaluating T Corr for different locations x gives a correlation map. Correlation-based methods offer optimal performance under Gaussian noise [START_REF] Kay | Fundamentals of statistical signal processing[END_REF]. It is well known, however, that the performance strongly deteriorates in the presence of large deviations (outliers), see [START_REF] Fitch | Fast robust correlation[END_REF] and figure 1.

T Corr (x) = max   m(x) t Γ -1 y m(x) t Γ -1 m(x) , 0   H1 ≷ H0 γ ,
Robust image matching is a widely studied topic. A common approach starts by detecting feature points, then computes local descriptors such as scale-invariant feature transforms (SIFTs) [START_REF] David G Lowe | Distinctive image features from scale-invariant keypoints[END_REF] before matching pairs of descriptors. Robust matching is achieved by restricting the matching to feature points, and by using well-chosen (i.e., invariant) descriptors [START_REF] Morel | ASIFT: A new framework for fully affine invariant image comparison[END_REF].

In order to detect patterns even at low signal-to-noise (SNR) regimes, a "pixel-based" approach is often better suited. The fast correlation method proposed in [START_REF] Fitch | Fast robust correlation[END_REF] for image registration replaces the quadratic penalty typical of maximum likelihood under a Gaussian assumption by a bounded penalty of the form ψ(∆) = 1 -cos(∆). The reason for this specific choice of penality is that it leads to a fast implementation based on FFTs.

In a signal processing context, locally optimal detectors (i.e., detectors statistically optimal for weak signals) have been developed for several noise distributions, see the seminal papers by Capon [START_REF] Capon | On the asmptotic efficiency of locally optimum detectors[END_REF] and by Miller and Thomas [START_REF] Miller | Detectors for discrete-time signals in non-Gaussian noise[END_REF]. These detectors have been applied recently to the detection of marks in image watermarking [START_REF] Bian | Locally optimal detection of image watermarks in the wavelet domain using Bessel K form distribution[END_REF], [START_REF] Briassouli | Locally optimum nonlinearities for DCT watermark detection[END_REF]. In this paper, we derive the general expression of locally optimal detectors for problem [START_REF] Kay | Fundamentals of statistical signal processing[END_REF], including the normalization of the test statistics required for the calibration of the threshold. This normalization is crucial to account for truncation at the image borders and prevent the introduction of bias for objects centered close to or beyond the image borders. It is also necessary when considering noise with a non-stationary dispersion. We show that the expression of the criterion offers a fast implementation based on pixelwise operations and discrete correlations. We give its closeform expression under Cauchy noise and prove its Gaussian approximation in the asymptotic regime. Finally, we illustrate its performance on a numerical experiment.

II. PROPOSED ROBUST DETECTOR

A. Motivation

Several tests have been developed in the framework of detection theory. The three most well known are the (generalized) likelihood ratio, Wald's statistic and Rao's score, see [START_REF] Rao | Score test: historical review and recent developments[END_REF], [START_REF] Erich | Testing Statistical Hypotheses[END_REF]. Robust detection criteria can readily be derived from these tests by modeling the distribution of errors p( ) with a heavy-tailed distribution. However, if the location x of the pattern is not known (at least approximately) beforehand, robust detection criteria are generally computationally very demanding because they require estimating the amplitude α during the exhaustive search for location x. Maximum likelihood estimation of α under non-Gaussian distributions (an Mestimator) typically requires resorting to an iterative procedure [START_REF] Peter | Robust statistics[END_REF]. Testing all pixel locations thus leads to prohibitive computation times when the pattern m has a large spatial extension.

A very interesting feature of Rao's score is that it does not require the estimation of the amplitude of the pattern. When the (unknown) amplitude of the pattern to detect is always nonnegative, Rao's score corresponds to the locally most powerful (LMP) detector, i.e., it is asymptotically optimal for weak signals (when α tends to zero). The LMP test is given by [START_REF] Kay | Fundamentals of statistical signal processing[END_REF]:

T LMP (y) = ∂ log p(y|α, H 1 ) ∂α α=0 • I -1/2 0 H1 ≷ H0 γ , (3)
where I 0 denotes Fisher information on parameter α, evaluated at α = 0:

I 0 = E - ∂ 2 log p(y|α, H 1 ) ∂α 2 α=0 . (4)
The LMP test does not require the estimation of the amplitude α, but only the evaluation of the slope of the log likelihood at α = 0 (i.e., under H 0 ) and the computation of the expectation of the curvature of the log likelihood for the Fisher information. In the case of Gaussian noise, the LMP detector leads to the same expression (2) as the generalized likelihood ratio. We investigate the use of this test with heavy-tailed distributions as a detection criterion robust to the presence of outliers.

The general expression of LMP detectors has been derived by Miller and Thomas [START_REF] Miller | Detectors for discrete-time signals in non-Gaussian noise[END_REF], without the normalization factor I 0 . It is worth to note that, as already mentioned in the introduction, the normalizing factor I 0 plays a crucial part in image processing. I 0 is a function of m(x) which generally depends on the position x. While the threshold γ can be calibrated once and for all, using Monte-Carlo simulations or asymptotic distribution of the test statistic, in the center of the image, i.e. where m is constant, a different calibration of γ will be required for each position of the pattern on the border of the image. Moreover, if the noise displays a non-stationary dispersion, the normalization I 0 also varies spatially.

B. Derivation of the detector

We define the scaled residuals t(α, x) by

∀i, t i (α, x) = y i -α m i (x) s i (5)
where we assume that the scaling factors s i > 0 can be chosen such that under H 0 , random variables t i (0, x) are independent and identically distributed (i.i.d.). We thus consider a slightly more general case than that of i.i.d. noise. By considering nonconstant scaling factors s i , one can account for noise with non-stationary dispersion. We further assume that the likelihood of the data can be written under the form

p(y|α, m(x)) ∝ i exp(-ϕ(t i )) , (6)
with ϕ twice differentiable.

To compute the first factor of the test (3), we begin by computing the partial derivative of the log-likelihood with respect to the amplitude α of the pattern:

∂ log p(y|α, H 1 ) ∂α = i m i (x) s i ϕ (t i ) , with ϕ = ∂ϕ ∂t (7)
which gives:

∂ log p(y|α, H 1 ) ∂α α=0 = i m i (x) s i ϕ (y i /s i ) . (8)
The second part requires the computation of Fisher information. We therefore compute the second order partial derivative:

∂ 2 log p(y|α, H 1 ) ∂α 2 = - i m 2 i (x) s 2 i ϕ (t i ) . (9)
Fisher information on parameter α, for α = 0, is given by:

I 0 (x) = i m 2 i (x) s 2 i E[ϕ (y i /s i )] . (10)
Since under our initial hypothesis t i (0) = y i /s i are i.i.d., the expectation is constant:

∀i, E[ϕ (y i /s i )] = ∞ -∞ ϕ (t) exp(-ϕ(t)) dt β , (11)
which gives the final expression of the LMP test:

T LMP (x) = τ (x) n i=1 m i (x) s i ϕ (y i /s i ) , (12)
with

τ (x) = β n i=1 m 2 i (x) s 2 i -1/2 . ( 13
)
We have written explicitly the dependency on the location x of the model in order to emphasize that the normalization coefficient τ (x) varies in the field of view, either because the model is shifted with respect to non-constant scaling factors s i , or because truncation and/or deformation of the model occur.

We further discuss the practical implication of the spatial variation of τ (x) in paragraph II-E.

C. Asymptotic normality under H 0

The model under consideration is not i.i.d. and consequently asymptotic normality of the LMP test statistic under H 0 cannot be assumed. The expression of the test obtained in equation ( 12) can be re-written under the form:

T LMP (x) = 1 n i=1 c 2 i (x) n i=1 c i (x) u i , (14)
with c i (x) = m i (x) s i and u i = ϕ (y i /s i )/ β . (15)
The expectation of random variables u i under H 0 is 0, while their variance is equal to 1. By application of lemma 11.3.3 from [START_REF] Erich | Testing Statistical Hypotheses[END_REF], a sufficient condition for the test statistic T LMP (x) to converge to a standard normal distribution under H 0 is that:

lim n→∞ n i=1 c 2 i (x) max i=1..n c 2 i (x) = ∞ (16)
which amounts to requiring that additional observations increase Fisher information so that it is unbounded.

A threshold can thus be set to reach a prescribed false alarm rate based on a normal approximation of the sample distribution of T LMP (x).

D. Derivation under Cauchy noise

Cauchy distribution has heavy tails that can account for the presence of strong outliers in the data. The negative loglikelihood for a standard Cauchy distribution is given by:

ϕ(t) = log(t 2 + 1) + log(π) , (17)
and its derivative is

ϕ (t) = 2t t 2 + 1 . ( 18
)
By application of equation ( 12), we obtain the expression of the robust detection test:

T LMP (x) = 2τ (x) n i=1 m i (x) s i y i /s i (y i /s i ) 2 + 1 , (19)
where the constant β in the normalization factor τ (x) can be computed in closed form [START_REF] Briassouli | Locally optimum nonlinearities for DCT watermark detection[END_REF]:

β = 1/2.
Note that since the derivative of the log-likelihood of Cauchy distribution is bounded by ±1, the random variables u i = ϕ (y i /s i )/ √ β are uniformly bounded (both under H 0 and H 1 ). Moreover, the following expressions for the mean and the variance of u i can be derived:

E[u i ] = 2 √ 2αs i m i (x) 4s 2 i + α 2 m 2 i (x) , Var[u i ] = 4s 2 i 4s 2 i + α 2 m 2 i (x) (20)
and the probability density function of u i is:

p(u i ) =      √ 2 2π 1 -4 ui 2 √ 2 -αmi(x) si 2 -1 2 if u i ∈] - √ 2, √ 2[ 0 otherwise. (21)
Based on (21), we can represent the distribution of the test statistic T LMP for several values of n. Figure 2 illustrates the convergence to the standard normal distribution under H 0 . Asymptotic normality of T LMP under H 0 can also be proved in the Cauchy noise case as a direct consequence of Lindeberg's theorem: every uniformly bounded sequence of mutually independent random variables obeys the central limit theorem, provided that the variance of the sum be unbounded (see [START_REF] Feller | An introduction to probability theory and its applications[END_REF]). Application of this result to the variables m i (x)u i /s i in (19) proves asymptotic normality of (19) when is also unbounded under H 1 . This proves that the sum of the m i (x)u i /s i , properly normalized, is also asymptotically normal under H 1 . This justifies a Gaussian approximation of the T LMP distribution under H 1 in the finite samples case. The asymptotic probability of detection and probability of false alarm can thus be computed based on the Gaussian cumulative distribution function.

E. Map of the robust-detection criterion

The expression obtained in equations ( 12) and ( 13) follows the general form of locally optimal (LO) detectors [START_REF] Capon | On the asmptotic efficiency of locally optimum detectors[END_REF], [START_REF] Miller | Detectors for discrete-time signals in non-Gaussian noise[END_REF]:

T LMP (x) = τ (x) • m(x) t g(y) , (22)
with g a non-linear but separable transformation applied on the vector of observed amplitudes y before correlating with the model m: ∀i, [g(y

)] i = 1 s i ϕ (y i /s i ) . ( 23
)
and in the special case of Cauchy noise, [g(y)] i = 2 y i /(y 2 i + s 2 i ). Input values y i that are small with respect to s i are left unchanged (up to a constant scaling factor), while larger values of y i are attenuated in order to reduce the influence of outliers. An attractive feature of the locally optimal tests for image processing is that they offer fast implementations based on pixel-wise operations and discrete correlations that can be computed using FFTs:

T LMP (x) = m (g(y) w) √ β (m 2 ) (w/s 2 ) ( 24
)
where w is a binary mask indicating the region where pixels are measured, s is a vector containing the scale factors s i , denotes the element-wise product, denotes the bidimensional correlation, and the divisions as well as the square root and squares are applied element-wise. Note that the nonlinear transform g is fast to compute since it is separable.

Figure 3 illustrates the importance of the normalization (i.e., denominator in equation ( 24)) when computing maps of the robust-detection criterion. Without proper normalization, the maximum of the criterion is always reached inside the field of view, which biases the detection.

F. Estimation of the s parameter in the case of Cauchy noise

The application of the robust criterion derived for Cauchy distribution requires the estimation of scaling factors s i . We consider in this paragraph the estimation of a constant parameter s. When the scaling factors s i are chosen constant and equal to s, the expression of the robust detector under Cauchy noise simplifies into:

T LMP (x) = 2 √ 2 m(x) 2 i m i (x) y i /s (y i /s) 2 + 1 . ( 25
)
In the context of the detection of weak signals, we suggest estimating s from the data by considering that noise is dominating, or that data where the pattern is not present are available. Maximum likelihood estimation of s is done by finding the unique root of the derivative of the log-likelihood, i.e., by solving the non-linear equation [13]:

i 2 s 2 s 2 + y 2 i = n ( 26
)
where n is the dimension of the observation vector y. Equation ( 26) can be solved efficiently by bisection, starting from the initial value for s: i |y i |/n, or by using Brent's method [START_REF] Richard | Algorithms for minimization without derivatives[END_REF]. This estimation algorithm is preferred to algorithms relying on order statistics, such as [START_REF] Kravchuk | Hodges-Lehmann scale estimator for Cauchy distribution[END_REF], for its computational efficiency.

III. COMPUTER SIMULATIONS

We illustrate the detection performance of the proposed criterion on a numerical simulation. The reference pattern is Gaussian (fig. 1(a)). The observations y are corrupted by Gaussian mixture noise: 95% of the pixels noise is Gaussian with an SNR of 1 and the remaining 5% is Gaussian with an SNR of 0.01 (outliers). We illustrate in figure 1(e) that the map of the detection criterion values defined in equation ( 24) reaches a maximum close to the true location of the pattern, while standard correlation, i.e., the detector defined in equation ( 2), gives a maximum at an incorrect location. As illustrated in figure 3, proper normalization is crucial to account for possible truncation of the pattern. If the normalization is missing, a bias in the localization appears.

The detection performance is further compared by computing the receiver operating characteristic (ROC) curve from 100000 random realizations in the same setting, figure 4. We compare the performance of several detectors: (i) the Neyman-Pearson idealized detector that uses the true distribution of noise (a mixture of two Gaussians) and the knowledge of the exact value of the amplitude α; (ii) the Generalized Likelihood Ratio (GLR) for a Cauchy distribution (parameter s is estimated from the data, even when α = 0, using the method described in paragraph II-F); (iii) the proposed robust detector, using the same estimation of Cauchy distribution parameter s; (iv) the proposed robust detector with the best possible tuning of s (i.e., which maximizes the area under the ROC curve); and (v) the usual correlation criterion given in equation ( 2). From the ROC curves, it can be observed that correlation performs poorly because of the presence of outliers. Robust criteria behave much better. The GLR and the proposed LMP criteria offer very similar performance, the latter being much easier to implement since it does not require to estimate α. Note that the distribution of noise in the numerical experiment (Gaussian mixture) differs from the hypothesized distribution in the derivation of the test, without significantly impacting the performance: with the best possible choice of s, the performance of our criterion almost reaches that of Neyman-Pearson's detector.

Beyond correct detection, it is generally necessary in applications to both detect and locate the pattern. The same noise model as in figure 4 is considered, only with a stronger SNR (95% of the pixels with an SNR of 8 and the remaining 5% with an SNR of 0.08). The pattern corresponds to the truncated Gaussian with a standard deviation equal to 40 pixels that is depicted in figure 4. The location of the pattern along the y-axis is estimated by finding (by exhaustive search) the location that maximizes each test statistic. The x-axis location is kept constant and equal to the correct location. The standard deviation of the location, computed over 1000 random trials, is reported in table I for each criterion. In all cases, we found the bias to be negligible compared to the standard deviation. The Neyman-Pearson criterion, based on the knowledge of the exact distribution of noise and of the true amplitude α of the pattern gives the best accuracy. Maximizing this criterion amounts to estimating the y location with the maximum likelihood estimator using the exact likelihood and knowledge of the amplitude α. Maximizing the normalized correlation, i.e., maximum likelihood estimation under the hypothesis of a stationary Gaussian noise, gives a very strong localization error (the error is multiplied by 29). Such a non-robust criterion is not adapted to detection and estimation problems in the presence of strong outliers. Maximizing the GLR under the hypothesis of a Cauchy distribution, with the parameter s of Cauchy distribution estimated from the data, slightly degrades the accuracy compared to the (idealized) Neyman-Pearson criterion. This estimation method corresponds to the joint maximum likelihood estimation of the location y and amplitude α under an hypothesized Cauchy distribution, i.e., an M-estimator. Finally, maximizing the proposed criterion (the LMP) gives a comparable degradation of the accuracy. The huge difference in practice is that maximization of the LMP statistic does not require to jointly estimate α. It can thus be performed very efficiently on a pixel grid by FFTs.

IV. CONCLUSION

We derived the general expression of a robust detection criterion based on the theory of locally optimal detectors. We have shown the importance of the normalization term for image processing applications and have given a fast algorithm to compute detection maps. Application of the criterion to Cauchy distribution leads to a simple expression. The scale parameter of Cauchy distribution can be estimated directly from the data and the obtained detector offers a significant improvement compared to classical correlation while retaining its computational efficiency.

Fig. 1 .

 1 Fig. 1. Detection of a pattern in the presence of strong outliers: correlation vs. robust correlation. (a) the reference pattern; (b) noiseless signal: measurements are available only inside the white square, the center of the actual pattern is marked with a red circle; (c) noisy data (SNR=1 for 95% of the pixels, and SNR=0.01 for the remaining 5%) shown using the same display range as image (b); (d) normalized correlation map T Corr (maximum marked with a star); (e) proposed robust-correlation map.

Fig. 2 .

 2 Fig. 2. Asymptotic normality of the test under Cauchy noise: probability distribution of T LMP for different values of n, under H 0 , with s i = 1 and m i = 1 for all i, and α = 1.

n i=1 c 2 i

 2 (x) = ∞.Note that when α > 0, if c i (x) is upper bounded, the variance (20) is lower bounded. Hence when n i=1 c 2 i (x) = ∞ the sum of variances of the variables m i (x)u i /s i in(19)

Fig. 3 .

 3 Fig. 3. Importance of proper normalization of the detection criterion: detection of a Gaussian pattern centered on the red circle, without normalization (second column) or with (last column) normalization.

Fig. 4 .

 4 Fig. 4. ROC curves for several detectors. The curves of the GLR and the proposed LMP test are almost perfectly superimposed.

TABLE I LOCALIZATION

 I ACCURACY OF SEVERAL CRITERIA

	Criterion	Localization accuracy (standard deviation in pixels)
	Neyman-Pearson	2.2
	Correlation	64.5
	Generalized Likelihood Ratio (Cauchy)	2.8
	Proposed criterion (Cauchy)	2.9

This work has been supported by a grant from the French CNRS (DETEC-TION, DEFI IMAG'In 2015)