Vladimir A Krylov
email: vladimir.krylov@unige.it

S. B. Serpico are with DITEN Gabriele Moser
email: gabriele.moser@unige.it

Fellow, IEEE Sebastiano B Serpico
email: sebastiano.serpico@unige.it

Fellow, IEEE Josiane Zerubia
email: josiane.zerubia@inria.fr

False Discovery Rate Approach to Unsupervised Image Change Detection

Keywords: Change detection, statistical hypothesis testing, false discovery rate, synthetic aperture radar images, dermatological images, still camera images. I

In this paper we address the problem of unsupervised change detection on two or more coregistered images of the same object or scene at several time instants. We propose a novel empirical-Bayesian approach that is based on a false discovery rate formulation for statistical inference on local patch-based samples. This alternative error metric allows to efficiently adjust the family-wise error rate in case of the considered large-scale testing problem. The designed change detector operates in an unsupervised manner under the assumption of the limited amount of changes in the analyzed imagery. The detection is based in the use of various statistical features, which enable the detector to address application-specific detection problems provided an appropriate ad hoc feature choice. In particular, we demonstrate the use of the rank-based statistics: Wilcoxon and Cramér-von Mises for image pairs, and multisample Levene statistic for short image sequences. The experiments with remotely sensed radar, dermatological, and still camera surveillance imagery demonstrate accurate performance and flexibility of the proposed method.

is immediate that in most cases a very particular application-specific kind of changes is sought for. In fact even the considered types of image sequences are very different: from static camera videos with tens of frames per second to two satellite images obtained with a time gap of several months and significant variation in acquisition conditions, like incidence angle and weather effects. When a specific change detection application is considered the first thing to determine is the type of significant changes that is to be detected. This is directly related to the considered class of disturbance factors, i.e., the reasons that bring to non-significant changes that should be ignored by the detector.

These might include acquisition noise, coregistration artefacts, illumination changes, etc. Naturally, there is an appreciable amount of research work addressing various change detection problems [START_REF] Radke | Image change detection algorithms: a systematic survey[END_REF], [START_REF] Hussain | Change detection from remotely sensed images: From pixel-based to object-based approaches[END_REF]: most popular techniques include the background subtraction [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF], [START_REF] Barnich | ViBe: A universal background subtraction algorithm for video sequences[END_REF], [START_REF] Haines | Background subtraction with Dirichlet process mixture models[END_REF], [START_REF] St-Charles | SuBSENSE: A universal change detection method with local adaptive sensitivity[END_REF], and local patch approaches based on pixel ordering [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF], [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF].

In this paper we focus on the classic change detection scenario when changes are evaluated from two or several coregistered 1 images obtained at different time instants. The central assumption is that the number of changed pixels should be small relative to the number of unchanged pixels, allowing the algorithm to identify areas exhibiting abnormal statistical behavior as changes. Such a scenario arises when limited changes are sought within limited quantities of data, e.g., in surveillance, remote sensing and medical image processing problems. To address this problem we develop a generic statistical approach, that takes pixel-level detection decisions based on statistical hypothesis tests performed on grouped image samples. Since the construction of representative temporal samples is unfeasible, we consider local patch samples obtained with a square sliding window. The resulting spatial samples are processed with the appropriate test statistics, that are referred hereafter as features. The large choice of such features allows to take into account various disturbance factors, and, thus, address particular change scenarios.

The developed empirical-Bayesian approach belongs to the class of large-scale inference problems, where thousands of individual inferences (tests) are performed simultaneously (one per every image pixel). This formulation leads to a difficult problem of adjusting the family-wise error (FWE) [START_REF] Conover | Practical nonparametric statistics[END_REF], [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF], which consists in designing a specific multistep procedure to control the amount of the so-called type I errors, or, in other words, false positives.

We address this problem with a false discovery rate (FDR) approach [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF], [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF], a fairly recent statistical method that replaces the classical significance testing by an error metric that controls the proportion of false positives among all detections. This metric is appropriate for testing a large number of hypotheses when the standard statistical testing protocol is too conservative. The FDR-based techniques bring to more powerful testing procedures [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF], [START_REF] Efron | Large-scale simultaneous hypothesis testing: The choice of a null hypothesis[END_REF], and more robust to the patch-based sample dependences [START_REF] Farcomeni | A review of modern multiple hypothesis testing, with particular attention to the false discovery proportion[END_REF].

The contribution of this paper is in the development of a novel statistically sound FDR-based approach to image change detection. This approach allows to adapt a wide variety of the statistical tests existing in the literature to the problem of patch-based image change detection. We demonstrate the ability of the developed approach to produce accurate application-specific detection results by an appropriate choice of the employed feature. Specifically, in case of image pairs we show the use of Wilcoxon and Cramér-von Mises statistics, that verify the equality of first order moments and the equality of the distributions, respectively. In case of multiple images we propose the use of the Levene multisample statistic, that tests the hypothesis of the equality of the sample variances. For validation purposes we consider (i) remotely sensed synthetic aperture radar (SAR) imagery, where the main disturbance factor is given by the speckle noise, (ii) dermatological images with approximate coregistration and varying illumination conditions, and (iii) still camera images with disturbance factors of varying illumination and shadows. A preliminary short version of this study was presented in [START_REF] Krylov | False discovery rate approach to image change detection[END_REF].

The paper is organized as follows: In Section II we briefly review relevant previous work on statistical image change detection. In Section III we formulate the FDR-approach and propose the approaches to null distribution and empirical density estimation. In Section IV we propose a novel feature-based change detector. In Section V we present several statistical features for two and multiple images scenarios. We report the experimental results on three types of image pairs in Section VI, and a semi-synthetic experiment on an image sequence in Section VII.

The paper is concluded in Section VIII.

II. RELATED WORK

For a comprehensive overview of change detection methods we refer the reader to surveys performed in [START_REF] Radke | Image change detection algorithms: a systematic survey[END_REF], [START_REF] Hussain | Change detection from remotely sensed images: From pixel-based to object-based approaches[END_REF], [START_REF] Haines | Background subtraction with Dirichlet process mixture models[END_REF], [START_REF] St-Charles | SuBSENSE: A universal change detection method with local adaptive sensitivity[END_REF], [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF]. The brief selection of methods presented in this section is intended to give some insight into the recent developments in statistical image change detection. Note that we leave aside such wide areas of change detection techniques as segmentation and object-based detection, machine learning methods [START_REF] Hussain | Change detection from remotely sensed images: From pixel-based to object-based approaches[END_REF].

We first present several relevant FDR techniques that have been employed in the image processing applications.

In [START_REF] Hopkins | A new source detection algorithm using the false-discovery rate[END_REF] the FDR-analysis is performed on single-band astronomical images to identify specific (source) areas based on the p-values estimated via an a priori assumption of standard normal distribution for pixel intensities. In [START_REF] Mchugh | Motion detection with false discovery rate control[END_REF] the processing is done on long temporal sequences (100 frames) from surveillance cameras with the null assumption of normality of temporal samples at each location. In [START_REF] Van | False discovery rate for wavelet-based statistical parametric mapping[END_REF] and [START_REF] Nguyen | False discovery rate control in magnetic resonance imaging studies via Markov random fields[END_REF] the authors apply FDR-techniques to the MRI/fMRI imagery to identify brain areas (voxels) reacting to visual stimuli and subject to the ageing effect, respectively. In the former, the FDR-testing is performed on the iteratively thresholded output of the wavelet-based statistical parametric mapping (with the underlying normalized Student's test) of the input. In [START_REF] Nguyen | False discovery rate control in magnetic resonance imaging studies via Markov random fields[END_REF] the authors applied FDR-inference to the z-test p-values to identify the departures from zero of the voxel means (under normality assumption) in the simulated study, and to the Wald-like statistic on the coefficient of the age factor in the linear regression model of the voxel intensities in case of the real data. The authors have also proposed to integrate a second thresholding rule based on an isotropic MRF term to promote spatial regularity of the detection. In [START_REF] Seo | Action recognition from one example[END_REF] FDR is used to find video volumes similar to a given video query by resorting to the trace-based distances with the corresponding p-values recovered from the Student's t-statistic.

We now review several general classes of change detection methods. The first wide class is given by background subtraction approaches [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF], [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF], [START_REF] Barnich | ViBe: A universal background subtraction algorithm for video sequences[END_REF], [START_REF] Haines | Background subtraction with Dirichlet process mixture models[END_REF], [START_REF] St-Charles | SuBSENSE: A universal change detection method with local adaptive sensitivity[END_REF], [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF]. These methods analyze the temporal consistency of image pixels in large image sequences. The second large class of methods rely on the spatial support by processing the local information at two time instants. Typically, the order statistics are used due to the fact that pixel orderings are invariant to monotonic transforms of the data, and also allow more generality with respect to different noise models.

The disturbance factors effect has been modelled as affine [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF] or more general order-preserving transformations [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF]. The central disturbance factors in these papers are the illumination changes and camera adjustments. The approach developed in [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF] tests the rank-consistency under an explicit noise model. Some contour-based change detectors have been developed [START_REF] Yokoyama | A contour-based moving object detection and tracking[END_REF], [START_REF] Yilmaz | Object tracking: A survey[END_REF]. Although capable of operating on only two images (of which one is considered as a background), these methods require estimates (variances) of the acquisition/camera noise, that are typically obtained from available stack of previously acquired frames.

In this paper we focus on a scenario when only several images are available for the change detection analysis. This represents a different problem for which many of the background-based methods, e.g., [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF], [START_REF] Haines | Background subtraction with Dirichlet process mixture models[END_REF], [START_REF] St-Charles | SuBSENSE: A universal change detection method with local adaptive sensitivity[END_REF], [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF], cannot be applied. Some methods, like [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF], [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], can be employed in two frames case given the noise variance estimates. For comparisons we will use methods proposed by Lanza et al. [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], that verifies local order consistency and formulates detection as a maximum likelihood nonparametric isotonic regression problem, and Xie et. al. [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF], that relies on an explicit statistical imaging model and likelihood ratio test of order consistency.

Considering SAR change detection problem addressed in the experimental validation, most of the existing methods build on explicit parametric models [START_REF] Moser | Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery[END_REF]. Typically, due to the multiplicative nature of the speckle, SAR change detection methods rely on image ratio modeling, e.g., generalized Gaussian [START_REF] Bazi | An unsupervised approach based on the generalized Gaussian model to automatic change detection in multitemporal SAR images[END_REF], or using lognormal-or Nakagamiratio distributions [START_REF] Moser | Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery[END_REF]. The detection maps are typically obtained via minimum-error thresholding algorithms [START_REF] Moser | Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery[END_REF], [START_REF] Bazi | An unsupervised approach based on the generalized Gaussian model to automatic change detection in multitemporal SAR images[END_REF]. Further methods applied to remote sensing data are based on fuzzy clustering [START_REF] Gong | Change detection in synthetic aperture radar images based on image fusion and fuzzy clustering[END_REF], likelihood ratio testing [START_REF] Lingg | A sequential framework for image change detection[END_REF],

stochastic geometry [START_REF] Benedek | Building development monitoring in multitemporal remotely sensed image pairs with stochastic birth-death dynamics[END_REF], and multi-layer Markovian models [START_REF] Benedek | Detection of object motion regions in aerial image pairs with a multi-layer Markovian model[END_REF]. For comparisons we will employ explicit statistical models developed by Moser et al. [START_REF] Moser | Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery[END_REF], that is based on lognormal-ratio modeling, and by Inglada and Mercier [START_REF] Inglada | A new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis[END_REF], that builds on thresholding of cumulant-based Kullback-Leibler distance between local patch-based samples and can be applied to arbitrary types of imagery.

III. FALSE DISCOVERY RATE APPROACH

When addressing the change detection problem as a set of pixel-based decisions on a stack of images, we face a so-called large-scale simultaneous hypothesis testing problem [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF], [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. This type of problems is characterized by thousands of simultaneous tests that are performed on data samples of the same nature. Operating on two coregistered images, we calculate a value of a selected statistical feature F at each pixel location. In this paper we employ the statistic features, for which distributions D F are known in a no-change scenario under the classical sample independence assumption. Working with local patch samples, we encounter the problem of dependence in the values of collected features, at least, at neighbouring locations. Furthermore, the correction to the individual test significance level in a multi-test is a sophisticated matter [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. Specifically, it is well known that when addressing a multiple inference problem, an adjustment needs to be made to preserve the significance level of the resulting test. The classical remedy given by the Bonferroni correction [START_REF] Conover | Practical nonparametric statistics[END_REF] consists in setting a new significance level of each individual test: If it is desired that the significance level of the group-test with n hypotheses should be α, then the corrected individual test's significance level should be set to α/n. Naturally, such a level of significance is overly conservative and impossible to obtain in practice (considering n 10 6 pixels), which inevitably results in zero-detection. Other common methods imposing FWE-control include Tukey-Cramer method, step-down procedure, Games-Howell method, see [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF] for more detail. Another common technique to obtain empirical p-values, given by non-parametric bootstrap [START_REF] Conover | Practical nonparametric statistics[END_REF], cannot be employed here since the total amount of tests (equal to the image size) is far greater than the size of the samples employed for each test (local patch size), see [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF].

Due to the above reasons the direct use of the p-values obtained feature distribution D F is unfeasible. To address the problem, we suggest to capitalize on the knowledge of a significant amount of feature observations-one per each pixel of the considered image stack-to learn more about the actual D F distribution. This becomes possible since we assume that the majority of at least 75% of the pixels have not undergone any significant change. To learn from this large set of non-changed data we adopt the local FDR approach [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF], [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF] to image change detection.

We begin by introducing the notion of z-scores. This allows the approach to be formulated in a unified way regardless of the employed form of D F and assume that the transformed feature values (z-scores) are normally distributed. Let D F denote the cumulative distribution function (CDF) of the feature F, then we can define a z-score of the i-th image pixel as:

z i = Φ -1 (D F (f i )), (1)
where Φ -1 is the inverse of a standard normal CDF, and f i are the observed values of the feature F (centered in i-th pixel). Note that if D F is non-parametric then the calculation of z-scores can be performed directly from the observed features as defined in (1) without resorting to any additional parameter estimation procedure. Indeed, most of the standard hypothesis testing procedures [START_REF] Conover | Practical nonparametric statistics[END_REF] lead to non-parametric (asymptotic) distributions for their statistics. In this paper we consider only non-parametric D F . It is immediate that in case of a normal distributed feature the z-scores are identical to the feature values.

We now formulate the testing problem in a Bayesian form. Suppose that the values of the feature statistic belong to two classes corresponding to unchanged (we will denote it the null hypothesis, H F 0 ) and changed (alternative H F 1 ), with prior probabilities p 0 and 1 -p 0 . We define the considered feature densities f 0 (z) and f 1 (z) depending on the class:

f 0 (z) density if unchanged (H F 0 ), with probability p 0 , f 1 (z) density if changed (H F 1 ), with probability 1 -p 0 .
The complete feature density is, then, a mixture f (z) = p 0 f 0 (z) + (1 -p 0 )f 1 (z). Thus, according to the Bayes theorem we obtain the a posteriori probability

Prob{unchanged | z} = p 0 f 0 (z)/f (z). (2)
The local false discovery rate is defined as:

lFDR(z) ≡ f 0 (z)/f (z). (3)
By ignoring the factor p 0 in (2), lFDR(z) can be considered as an upper bound for the probability of observing the unchanged class for an observation with a given value of z-score. The standard FDR(z) is a Bayesian-based formulation of a frequentist FDR quantity, and is expressed as an expectation of lFDR(z) conditioned on Z z [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. The local FDR is the probability that a particular test gives rise to a false positive. The advantage of lFDR(z)-control is that the independence is not required of the z-scores, all that is needed is a reasonable estimate of their marginal distribution.

In accordance with the considered change detection problem (assumption of relatively small amount of changes), we can assume that p 0 is sufficiently close to 1 (e.g., p 0 = 0.8 or p 0 = 0.9), so that lFDR(z) is a slightly conservatively biased estimator [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF] of (2). The FDR-based decisions in the following are taken by performing the thresholding and reporting z-scores with lFDR(z) γ as changed. This means that we label as changed the pixel locations that are the centres of local patches employed to produce the z-scores that are identified by the inequality above.

At this point it is important to underline the difference between a FDR-approach and more classical FWEprocedures. Both are designed to control the rejection of a true H F 0 hypothesis probability, however, whereas the FWE-methods control this quantity directly, the FDR-methods control the proportion of false positives among all rejected null hypotheses (discoveries). The advantage of the FDR-control is that it leads to fewer type II errors (false negatives) than controlling the FWE at the same level [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. We stress furthermore that the proposed FDR-based procedure is different from the conventional thresholding (possibly with different thresholds applied at two tails) of the histogram of z-scores which constitutes the standard statistical test. Specifically, the thresholding of the FDRfunction defined in (3) allows to consider pixels corresponding to some local z-value peak as changed, and at the same time leave less frequent z-values around in the unchanged class. From the statistical point of view that would mean that the z-values in the peak have been observed too often to consider them as following the null distribution.

Thus, this particular z-value and its neighbourhood are considered characteristic of one of the change scenarios.

A. Null distribution estimation

We proceed to the empirical DF 0 null distribution estimation based on the observed z-scores. This strategy allows accommodating the theoretical D F 0 distribution to the image patch samples, thus, addressing the dependence problems.

The most sensitive part in the considered methodological approach falls on the choice of the empirical null distribution f 0 (z) for the unchanged z-scores. If chosen perfectly, i.e., so that all the unchanged scores have f 0 (z) > 0 and f 1 (z) = 0 and the changed z-scores satisfy f 0 (z) = 0 and f 1 (z) > 0, then p0 is an unbiased estimate of the true percentage of the unchanged observations [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. Unfortunately such a choice is impossible in practice and, therefore, an empirical procedure needs to be established to estimate the null distribution f 0 (z).

As has been observed in typical applicative problems [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF], under the null hypothesis, the z-scores demonstrate a pronounced normal behavior even in cases, when substantial dependences are present in the data. Nevertheless, the classical standard N (0, 1) does not reflect accurately the statistical behaviour of z-scores at null even in conventionally dependence-free large scale inference problems [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. Therefore, in this paper we propose to consider the normal distribution family for null distribution, but allow the distribution parameters (mean and variance) to deviate from the standard N (0, 1) configuration in order to better reflect the dependence structure present in the data.

To arrive at the estimate of the null distribution that is affected by the changed pixels as little as possible we perform its estimation on the z-scores that lie around its distribution center and assume that the changed pixels report z-scores which lie far from the observed distribution center. In other words, the zero assumption is that the non-null density f 1 (z) is zero on a continuous interval containing the origin z * = 0, which renders the mixture model identifiable from the central histogram counts. This assumption holds if the chosen statistical feature is appropriate for the considered change detection problem, i.e., if it corresponds well to the significant changes. We employ the central matching approach [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF] by assuming that near the central score z * the representation f (z) = p 0 f 0 (z) holds, and that f 0 (z) follows the normal distribution N (δ 0 , σ 2 0 ). Under independence assumption z * should be set to zero, in practice we set it equal to the median of z-values. By taking the logarithm of f (z), we obtain

log f (z) = - 1 2σ 2 0 z 2 + δ 0 σ 2 0 z + log p 0 - δ 2 0 2σ 2 0 - 1 2 log(2πσ 2 0 ) .
The central matching strategy consists in estimating f 0 (z) by assuming that log f (z) is quadratic near z * , and, therefore, can be presented as:

log f (z) = β 2 z 2 + β 1 z + β 0 . (4)
We can thus obtain the estimates

( δ0 , σ2 0 ) = - β 1 2β 2 , - 1 2β 2 . (5)
Using ( 5) we arrive at an estimate of the empirical null f 0 (z) as N ( δ0 , σ2 0 ). A critical issue is the choice of the z-scores near z * , that are used to estimate (β 1 , β 2 ), and correspond only to unchanged observations. In this paper, we employ a conservative strategy by taking ∆ = 50% of the data for the fitting of (4). To define the interval of the feature histogram counts that are employed for f 0 (z) estimation, we adopt the following approach. We start by uniformly partitioning the range of z-values into K bins (intervals):

Z = K k=1 Z k . (6)
The bins are of equal widths, the number of z-scores in every bin is y k = #{z i ∈ Z k }, and x k are the centerpoints of Z k .

In order to construct the set of bins Z null used for the null estimation we first initialize it with the bin that contains the median z * . We then proceed by iteratively adding the adjacent bin (to the left or right from Z null ) that has the higher y k to the current set. This is iterated until the portion of scores inside the selected bins reaches 50% of the total:

Z k ∈Z null y k 0.5 K k=1 y k .
Then, the set of points (x k , y k ) characterizing the bins in Z null is fitted by the right hand side of (4) using a least square error approach [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. Finally, we get the estimates ( δ0 , σ2 0 ) as defined in [START_REF] Haines | Background subtraction with Dirichlet process mixture models[END_REF]. A typical example of such an empirical fit with f (z), f 0 (z) and N (0, 1) is presented in Fig. 1.

Note that the f 0 (z) estimates based on more than ∆ = 50% of observations might be more accurate, at the expense of potentially getting biased from the inclusion of changed z-scores, and, therefore, should not be considered as a default option. According to this strategy we define the above mentioned fraction of the changed pixels as limited by 25% of all pixels. This typically ensures that the central part of the appropriate feature histogram is well separated from the changed z-scores in the tail parts, and, thus, the change detection can be performed in an unsupervised manner. Furthermore, this assumption ensures the consistency of the obtained FDR-estimates under some technical assumptions. In particular, the Gaussianity of the central z-scores provides consistency of the least square estimate and the continuity of the exp-transform preserves this property. The consequent demonstration of the consistency for the FDR-estimates follows the proof of Theorem 1 in [START_REF] Nguyen | False discovery rate control in magnetic resonance imaging studies via Markov random fields[END_REF].

Alternative estimates can be obtained by maximizing the pseudo-likelihood (PL) of f (z), see in [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF], [START_REF] Nguyen | False discovery rate control in magnetic resonance imaging studies via Markov random fields[END_REF]. This procedure allows to obtain robust estimates and exploits the whole set of z-scores. Computationally, the Expectation-Maximization-like procedure necessary for the PL evaluation is more demanding. Experimentally, we have observed similar performance of central matching and PL, with the exception of datasets featuring more changes (up to 25%), where PL reported slightly more biased estimates of p 0 . For these two reasons we present the central matching procedure in this paper.

B. Empirical density estimation

To evaluate the FDR as defined in (3) we require a smooth estimate of the density f (z) for the observed values f i .

One way to address this problem is by resorting to classical interpolation tools such as natural splines [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. However, such tools do not take into account the normalization and non-negativity restrictions imposed on probability density functions. In this work, instead, we employ the approximation by a J-parametric exponential density function:

f (z) = exp J j=0 α j z j , (7)
where the constant α 0 is required to ensure the normalization. The choice J = 2 makes (7) into the normal family.

If J is large the estimate becomes almost non-parametric and tries to fit every z-value. A choice of J * = 7 in [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF] allows to achieve smooth estimates and does not impose too many restrictions on the density function shape.

To estimate the model parameters it is recommended to use the Lindsey's method that allows to obtain maximum likelihood estimates using a textbook Poisson regression. This involves the binning of the data as in [START_REF] St-Charles | SuBSENSE: A universal change detection method with local adaptive sensitivity[END_REF], and builds on the multinomial representation of the binned exponential statistics in [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF]. As such, the problem is conveniently reduced to a standard Poisson generalized linear model as described in [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF].

Algorithm 1: FDR-based change detector 1 calculate feature f(i, j) for all pixels (i, j);

2 transform f(i, j) into scores z(i, j) as in (1);

3 estimate the null density f 0 (z), see Sec. III-A;

4 approximate the empirical density f (z), see Sec. III-B;

5 calculate FDR(z) as in (3);

6 report as changed the pixels with FDR(z) < γ.

IV. IMAGE CHANGE DETECTION

In this section we formulate a novel feature-based change detector that operates on two or more images. We first assume that the significant changes are present in a fraction of the images' pixels that is at most 25% of their total size. We also recall that to employ the FDR-approach we can only use the statistical features with a known distribution under a certain null no-change hypothesis. Typically, any statistic that is used to construct a statistical hypothesis test can be employed as a feature in the designed change detector.

In order to accommodate the sample-based approach to construct the pixel-based detection map we employ the local patches. Specifically, for a pixel at location (i, j), we consider the sample obtained from the values of the surrounding pixels in a S × S local window centered at (i, j). It is immediate that any sample constructed by taking pixel values in a local-window of a natural image cannot be considered independent in itself, because of the contextual dependences associated with the image (spatial dependence or intra-sample dependence). Moreover, we consider the change detection problem, in which we assume that the amount of significant changes cannot be large.

As such, the two patch samples centered at pixel (i, j) on two coregistered images X and Y cannot be considered independent, i.e., for all i = 1, . . . , N observations at the same positions x i and y i are likely to be dependent As has been analysed in [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF], the FDR-formulation is robust to a certain level of individual tests' correlation for multiple inference. Considering the intra-sample dependence, we will assume that its level is moderate, so that the departures from statistic's distribution under the null hypothesis do not bring to appreciable non-normal behavior for the z-scores. Note that it is a common assumption for various applied problems, and that the normal statistical behavior of z-scores is commonly observed [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. Also note that the majority of existing change detection approaches ignore these dependence issues by working under various independence assumptions, see [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF], [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF], [START_REF] Moser | Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery[END_REF], [START_REF] Inglada | A new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis[END_REF].

The outline of the FDR-based change detector is presented in Algorithm 1. The choice of γ for the FDR-based procedure is similar to that of the significance level in classical inference. Finally, the local patch size S should be chosen manually with respect to the scale at which the significant changes manifest, which constitutes the same strategy as in [START_REF] Ohta | A statistical approach to background subtraction for surveillance systems[END_REF], [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF]. The computational complexity of the proposed algorithm is O(S 2 A log S): linear with respect to the number of the input pixels A and S 2 log S 2 -to the patch size S (due to the use of order statistics that rely on sample ordering with merge / block sorting). In this respect the algorithm is similar to [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF] and [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF], and faster than that in [START_REF] Inglada | A new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis[END_REF].

V. STATISTICAL FEATURES

The developed change detector constitutes a general approach that can be adopted for specific applications by an appropriate choice of statistical features. These should be the multisample statistics with a known distribution under the no-change hypothesis.

In this paper we focus on the use of non-parametric statistics F that are aimed at change detection from multitemporal images and allow a sufficient level of generality without restrictive distribution assumptions. Such generality is important when a single change detector is employed to reveal changes on imagery affected by several distinct disturbance factors. In particular, encouraged by the positive results obtained in [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF] we employ the order (rank) based statistics in the role of features. The rank transform is widely used in statistics to convert parametric tests into non-parametric by replacing the observations by their ranks [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF], [START_REF] Conover | Practical nonparametric statistics[END_REF]. The strict specific distribution assumptions in such cases are replaced by a very mild assumption of the similar distribution form.

Note, however, that the proposed technique allows the use of parametric statistics as well when there is a strong evidence that the analysed data follow a particular statistical model.

We consider the use of several features to address disturbance factors given by coregistration artefacts, acquisition noise and illumination variation. Namely, we use (i) the two-sample Wilcoxon paired statistic [START_REF] Conover | Practical nonparametric statistics[END_REF], that is decisive with respect to the first order sample changes, (ii) the two-sample Cramér-von Mises statistic [START_REF] Conover | Practical nonparametric statistics[END_REF] for arbitrary type of distribution changes, and (iii) the Levene statistic on the rank-transformed multisample data [START_REF] Conover | Practical nonparametric statistics[END_REF], [START_REF] Conover | Rank transformations as a bridge between parametric and nonparametric statistics[END_REF] that detects sample variance changes. These features are common in various applications thanks to their simplicity and robustness [START_REF] Conover | Practical nonparametric statistics[END_REF].

A. Two-sample Wilcoxon statistic

Given two independent random variables X and Y, the paired Wilcoxon statistic [START_REF] Conover | Practical nonparametric statistics[END_REF] verifies the null hypothesis

{H W 0 : P (X < Y) = 1/2}
, which characterizes the comparative central symmetry of two observed populations. The information given by this statistical feature characterizes the first order moment changes [START_REF] Conover | Practical nonparametric statistics[END_REF].

The two-sample paired Wilcoxon statistic is based on the observed absolute differences L i = |X i -Y i |, where X i and Y i are the observed samples of X and Y, respectively, i = 1, . . . , N . In the employed sliding window-approach N = S 2 , where S is the local patch size. The values of L i are transformed into the order statistic

L (1) < • • • < L (N ) .
Then the ranks of L (i) that correspond exclusively to the positive differences X i -Y i are summed to get W + . The standardized Wilcoxon statistic is defined as:

W = W + - N (N + 1) 4
N (N + 1)(2N + 1) 24 -0.5 (8)
This statistic converges in distribution to the standard normal N (0, 1) [START_REF] Conover | Practical nonparametric statistics[END_REF] when N → ∞ under H W 0 with the independence assumptions. Thus, for large samples, the critical values for the Wilcoxon test can be drawn from the normal approximation, that is classically employed when N > 20 [START_REF] Conover | Practical nonparametric statistics[END_REF]. The rejection region of the test is twosided, i.e., H W 0 is rejected when |W | > q, where q is a relevant quantile of the standard normal distribution. Since W ∼ N (0, 1), the statistic values f = W can be considered directly as z-scores, following the definition in [START_REF] Radke | Image change detection algorithms: a systematic survey[END_REF].

The paired Wilcoxon feature formulation (considering of differences |X i -Y i | instead of X i and Y i values) allows to compensate the inter-sample (temporal) dependence problem. The considered hypotheses H W 0 and H W 1 render the feature sensitive to the illumination changes (linear shifts), which, however, can be an important source of significant changes, e.g., when the central disturbance factor is given by zero-mean additive noise.

B. Two-sample Cramér-von Mises statistic

The two-sample Cramér-von Mises (CvM) statistic [START_REF] Conover | Practical nonparametric statistics[END_REF] is designed to verify a general null hypothesis: {H CvM Thus, compared to the Wilcoxon feature, the CvM feature can detect more complex changes like fluctuations in higher order moments. Note that similar results can be obtained with the other general two-sample goodness-of-fit features, e.g. Kolmogorov-Smirnov and Anderson-Darling statistics [START_REF] Conover | Practical nonparametric statistics[END_REF].

To calculate the two-sample CvM statistic, the observed samples X i and Y j are initially ordered into X (i) :

X (1) < • • • < X (N ) and Y (j) : Y (1) < • • • < Y (N )
. The ordered samples are concatenated into the sample L of size 2N , which is then ordered as well to obtain L (1) < • • • < L (2N ) . Then,

T = N i=1 (r(i) -i) 2 + (s(i) -i) 2 2N 2 - 4N 2 -1 12N (9)
defines the standardized CvM statistic, where r(i), s(i), i = 1, . . . , N , are the observations' ranks in the ordered samples: L (r(i)) = X (i) , and L (s(i)) = Y (i) . T follows asymptotically the Cramér-von Mises distribution [START_REF] Conover | Practical nonparametric statistics[END_REF], [START_REF] Xiao | A C++ program for the Cramér-von Mises two-sample test[END_REF].

However, the expression for the asymptotic CDF requires numerical approximations (since it involves the modified Bessel function of the third kind [START_REF]Handbook of Mathematical Functions[END_REF]), and the convergence of the distribution (as N increases) is relatively slow, especially so in the tail areas [START_REF] Xiao | A C++ program for the Cramér-von Mises two-sample test[END_REF]. For these reasons we will employ the exact CDFs for the sample size N as evaluated in [START_REF] Xiao | A C++ program for the Cramér-von Mises two-sample test[END_REF]. The rejection region of the CvM test is right sided. The z-scores can be estimated as defined in (1), and f = T .

Along with CvM, we propose the use of a modified CvM (mCvM) feature. It is constructed by preliminary subtracting sample medians from both samples before constructing T , see eq. ( 9). This allows to tune the statistic to fluctuations of higher order moments, ignoring the first moment changes, for example, constant illumination shifts. The distribution and critical region are the same as for CvM. Various further statistical features can be employed for detecting changes of higher order moments [START_REF] Conover | Practical nonparametric statistics[END_REF].

The non-paired CvM and mCvM construction allows to address the coregistration disturbance factor, since T is insensitive to the displacements of observations within large patches. The mCvM feature can be employed to address the illumination changes disturbance factor.

C. Multisample Levene statistic

In the case of change detection from a sequence of K images, consider K random variables X 1 , . . . , X K modeling the intensities in the same patch in K images. These variables are observed by samples of size N each:

(X k1 , . . . , X kN ) ∼ X k , k = 1, . . . , K. In this scenario there are two conceptual ways to validate the hypothesis that a certain characteristic of the samples is unchanged. The first consists in running an omnibus test that verifies the equality of all K sample characteristics. The term omnibus relates to the test formulation admitting multiple alternatives. Specifically, the null assumption of equality can be violated by any pair of samples which, when formulated explicitely, is equivalent to multiple distinct alternative hypotheses. The omnibus testing procedure is relatively fast, but in case when the null hypothesis is rejected the omnibus test does not specify in which pair(s) the changes occur. The second way is more computationally intensive and requires running all the pairwise tests.

These two options will be explored in this section in case of the Levene statistic evaluating the hypothesis of equal population variances. Note, that this statistic can be readily modified to test the equivalence of the sample means [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF], and allows further adjustments to target higher order statistical characteristics [START_REF] Rayner | Extended anova and rank transform procedures[END_REF].

1) Omnibus Levene statistic:

The null hypothesis is that the population variances σ 2 k are not changing: {H Levene 0 :

σ 2 1 = . . . = σ 2 K }.
The omnibus alternative stipulates that there exist two distinct populations X u and X v such that σ 2 u = σ 2 v . Note, that this test can also be employed in the two-sample scenario. The Levene test [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF], [START_REF] Nordstokke | A new nonparametric Levene test for equal variances[END_REF] is formulated as a repeated measures analysis of variance (ANOVA) test on the absolute deviation of the scores X kn from their corresponding sample medians Xk : Levene ≡ ANOVA |X kn -Xk | .

The classical formulation of the Levene test assumes the normality of the variables under scrutiny. One solution to relax this tight condition and tune the test to the non-parametric scenario exlored in this paper is provided by the concept of rank transformation [START_REF] Conover | Practical nonparametric statistics[END_REF]. This approach consists in performing the test on ranks rather than on the initial observations themselves, i.e.:

ANOVA |R kn -XR k | , (10)
with R kn denoting the rank of X kn in the pooled sample constructed by merging obervations from X 1 , . . . , X K with KN observations. Rk is the median rank in the k-th sample. Whenever ties occur, the average ranks are used.

Such a strategy constitutes a statistically valid and efficient testing procedure for multisample problems [START_REF] Nordstokke | A new nonparametric Levene test for equal variances[END_REF], [START_REF] Zimmerman | Relative power of the wilcoxon test, the friedman test, and repeated-measures anova on ranks[END_REF]. 2 As an omnibus test, [START_REF] Conover | Practical nonparametric statistics[END_REF] has slightly less power than single-alternative tests and is more sensitive to outliers [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF].

Therefore, to reduce the impact of noise and outliers on the detection results, we will take measures to render the test more robust by clipping [ωN ] observations in each of the input rank samples. Here ω is a predefined clipping percentage which is typically set to small values like 0.05, 0.1, see in [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF], and [x] denotes the integer part of x.

The Levene statistic for the dependent samples is given by

F = Q c /(K -1) Q e / (N -[ωN ] -1)(K -1) , ( 11
)
2 In particular, compared to the textbook non-parametric Friedman test [START_REF] Conover | Practical nonparametric statistics[END_REF], which is the multisample counterpart of the sign test, Levene test has the following advantages: (i) higher values of the asymptotic relative power [START_REF] Zimmerman | Relative power of the wilcoxon test, the friedman test, and repeated-measures anova on ranks[END_REF], and (ii) more sensitive p-values due to the ranking procedure performed on KN observations, and, thus, (KN )! possible combinations, which is substantially larger number than (N )! K combinations available for the Friedman test. where the estimation of quantities Q c and Q e from the ranks R kn is presented in Appendix A.

The critical values for the statistic F are taken from the Fisher CDF F v1,v2 with degrees of freedom v 1 and v 2 whose evaluation is also reported in Appendix A. The rejection region is right sided. The CDF of the Fisher distribution is expressed as an incomplete beta function, whose values can be easily obtained numerically [START_REF]Handbook of Mathematical Functions[END_REF] in order to produce the z-scores.

2) Pairwise Levene testing: Once the Levene test rejects the hypothesis of constant variance, one could ask a question which pairs of samples (images) in the sequence demonstrate meaningful changes and which not. To answer this question a panel of pairwise comparisons is required: if the omnibus test considered K populations, then K(K -1)/2 pairwise tests need to be performed. This procedure is commonly referred to as post-hoc testing [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF] because it relies on the result of the first test. To adjust the FWE in such multiple testing procedure we suggest to employ the Benjamini-Hochberg procedure [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF], which has more power than the classical Bonferroni correction.

This brings to a frequentist FDR inference procedure presented in Algorithm 2. The aim here is to ensure than the FDR is less or equal than γ. The rejected hypotheses correspond to pairs of images that present statistically meaningful differences in the variance parameter. Note, that in order to obtain the FWE-correction, an alternative

VI. EXPERIMENTS WITH IMAGE PAIRS

In this section we validate the proposed FDR-based approach on remotely sensed SAR, dermatological and still camera image pairs. The threshold for the FDR, which is the only conceptual parameter of the model and needs to be selected a priori, is set to γ = 0.1, see discussion in [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF]. However, examples of receiver operator characteristic (ROC) curves of the proposed method as a function of the threshold γ are also reported. The central matching is performed on ∆ = 50% of the data z-scores, and the partitioning with K = 75 provides good smoothness of the estimation. For the benchmark methods [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF] and [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF] the variance levels have been set in supervised manner in remote sensing and dermatological applications, estimated from 100 frames for the still camera application, and the thresholds (for Fig. 6) were set manually by choosing the best test performance. The experiments were performed on a C++ implementation with CPU-parallelization on an Intel Core-i7 2GHz, 8Gb RAM, Win7 system.

A. Remotely sensed images

We first present the experiments with speckled remotely sensed SAR imagery. Speckle is a noise-like phenomenon that appears as a result of mutual interference of wavefronts. It is characteristic of coherent imaging, in particular, in radar, sonar, lidar imagery [START_REF] Oliver | Understanding Synthetic Aperture Radar Images[END_REF]. Here we consider a flood monitoring sequence of two radar images in an agricultural area near Pavia, Italy, see Fig. 2, acquired by the shuttle-born XSAR sensor on (a) April 16 and (b)

April 18, 1994 (approx. 15 m res., four-looks). The performance was evaluated using a manually prepared ground truth map, the main changed areas are highlighted in Fig. 2(b). It is immediate, that the results obtained by the proposed method, with (e) CvM and (f) Wilcoxon features, are accurate and that they are less noisy than that of (d) the cumulant-based Kullback-Leibler (KL) thresholding [START_REF] Inglada | A new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis[END_REF] and (c) SAR-specific lognormal distribution-based generalized Kittler-Illingworth thresholding (Ln-GKIT) [START_REF] Moser | Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery[END_REF] methods. The numerical accuracies and computation times are reported in Table I, and the ROC curves presented in Fig. 3(a). The FDR is the expectation of the false discovery proportion (FDP) quantity presented in Table I. A ROC curve consists of only one point for Ln-GKIT method, since there is no input parameter; and for the proposed FDR-based methods, the γ parameter has been considered as a variable in order to construct its ROC curve. We observe that the employed features report more accurate results and are stable to the speckle disturbance factor with a slight oversmoothing on borders as a consequence of the employed sliding window approach.

Another remotely sensed SAR image pair, obtained by the space-born COSMO-SkyMed (CSK) sensor, depicts the port area of Port-au-Prince, Haiti, (a) before (April 2009) and (b) after (January 2010) the earthquake of January 2010 (StripMap acquisition mode, 2.5 m resolution, single-look), see Fig. 4. The purpose of this experiment is to compare visually the quality of the detection in a complex urban context. We report the benchmark results by (c) KL-based detection with thresholding at 2%, (d) Ln-GKIT [START_REF] Moser | Generalized minimum-error thresholding for unsupervised change detection from SAR amplitude imagery[END_REF], and (e) Lanza et al. method [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF]. In (f) we demonstrate the results of the classical inference procedure (α = 0.01 for individual tests) based on p-values of the Wilcoxon test, ignoring the intra-sample dependence. We then present the results obtained by the proposed FDR-based detector with: (g) the Wilcoxon feature, and (h) the mCvM feature. In the absence of ground-truth data we qualitatively observe that the results obtained with these different features are substantially different. The mCvM result is a more conservative version (with less detection) of the SAR-specific approach result reported in (d), both highlighting the changes in urban areas. Meanwhile the Wilcoxon result in (g) is more similar to general change detector result reported in (c). The visual examination suggests less noisy detection by the proposed FDR-based method as compared to methods in (c)-(e).

B. Dermatological images

We next validate the designed change detector in the application to hyper-pigmentation clinical treatment progress detection, see Fig. 5. Initial acquisitions were provided by one of the world leading companies in dermatology. The images were obtained by an optical sensor and evaluated by medical experts to identify the degree of significant change (one number per image pair). The images were coregistered and the zones of interest were defined (inside the inner contours): Solely the pixels inside these zones are tested. To compensate for the strong illumination changes the statistically meaningful thresholding procedure whereas the benchmark methods [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF] only produce profiles. Furthermore, for this specific quantitative application the control of false detection given by FDR is a strong advantage.

C. Still camera images

Here we present the application of the developed method to a classical still camera surveillance application [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF] on the images from the changedetection.net database [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF], see Fig. 6. The detection results are compared with that of the benchmark order-based methods proposed by Xie et al. [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF] and Lanza et al. [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF]. The mCvM feature (a)-(c), with γ = 0.05, gives a contouring detection, similar to the benchmark methods. In case of the mCvM feature, this effect is due to its construction that is sensitive to all changes (since Cramér-von-Mises test is derived to validate the general goodness-of-fit hypothesis), except for that of the first moment, which is subtracted from the individual samples. Consequently, when one homogeneous area is changed into another characterized by a distinct mean value but similar shape of the distribution (variance, and other higher order characteristics), the change is detected only along the border, which is the effect of patch-based formulation. Furthermore, the greater the considered patch size, the wider border is detected along the areas with an abrupt mean-change.

The image pair (c) is obtained by manually varying the brightness of the images in order to demonstrate the robustness of the mCvM feature to the constant intensity shifts. The above mentioned bordering effect is observed both on significant changes (moving people), and along the illumination shift border. This effect is, however, mitigated by the (partial) detection of the inner parts of the moving shapes, and, as such, the bordering is less pronounced than what is obtained with both [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF] and [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF]. Interestingly, the method in [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF] did not produce false positives along the intensity shift, instead failed to detect any significant changes inside this area. This might be due to the single level of variance set for the whole frame that resulted in the reduced sensitivity in the lower intensity area. The Wilcoxon (d) and CvM (e) features detect the inner-parts of the changed objects as well, yet

VII. EXPERIMENTS WITH MULTIPLE IMAGES

We present semi-synthetic experiments to validate the performance of the change detector in a three-image scenario. To this end we have manually constructed semi-synthetic 1000-by-1000 pixels SAR images from 6 zones, see Fig. 7: unchanged areas A and B that represent the vast majority of the pixels to ensure the fulfilment of the We present detection results in Fig. 7(e)-(h). Throughout this experiment we have employed sample clipping with ω = 0.05, and a 9-by-9 local window. First, in (e) we report an FDR map of pairwise detection of changes between images 2 and 3, and in (f) a combined three-band FDR map of all pairwise comparisons in inverse colour (yellow corresponds to image pair 1-2, sky blue -2-3, magenta -1-3, black -to changes in all three pairs). Note that here, as well as in (h), we perform soft FDR thresholding (i.e., do not replace the values below the FDR threshold) in order to demonstrate the contrast provided by the Levene statistic. In (g) we report the pre-FDR thresholding p-values collected from the null estimated using [START_REF] Haines | Background subtraction with Dirichlet process mixture models[END_REF], and in (h) the FDR-based detection using the omnibus Levene statistic.

The results in (e), (f) and (h) are post-processed by applying morphological reconstruction by opening [START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF]. It is employed with the smallest possible 3-by-3 filter window and serves to remove the isolated salt-and-pepper-like detected change pixels.

The analysis of the results demonstrates that the pairwise tests are more sensitive to changes, and produce results non-trivially different from those reported by the omnibus statistic. Indeed, the dark pixels in zone A in (e) correspond to small changes that are not distinguished by the omnibus statistic, see in (g). This omnibus decision is supported by the relative weakness of the change in the corresponding area identified by the other two pairwise tests. The same effect can be seen even better in zone IV in (f). Indeed, we observe a well visible change pattern in image pair 1-3, that is due to a strong change of distribution parameters used to generate the corresponding patches.

This detection is erroneous, since the distribution variance is in fact the same. The reason for that is the sample size which needs to be larger in order to adequately adjust to strong variation in the distributions. The omnibus test did not suffer from the same mistake due to a generally higher level of tolerance to pairwise fluctuations.

The performance in zones B and I-III is correct and rather comparable. Note, that here we have not relied on the post-hoc testing procedure presented in Algorithm 2, instead we have performed three pairwise and one omnibus Levene FDR detection procedures independently in order to evaluate their comparative performance.

Finally, we have empirically validated the stipulated above assumption of normality of z-scores with the dependencies present in the employed samples. We have observed Gaussian-like shape of the distribution even in the semi-synthetic experiment, which is less obvious since the images are comprised of areas of different nature-real and simulated SAR. Specifically, for the omnibus statistic the empirical null is estimated as N (1.61, 2.14) with Sk = 0.55 and Ku = 3.44, where Sk is the skewness and Ku -the kurtosis of the data [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF]. The latter two are typical characteristics of the normal distributions, that are exactly equal to 0 and 3, respectively, for perfectly Gaussian populations. In the pairwise analysis, the pair 2-3 demonstrated N (1.24, 2.55), Sk = 0.41, Ku = 2.85, and essentially similar values for the other two pairs. The values of skewness and kurtosis are not too far from the Gaussian values, and typically improve with larger sample sizes (local window). Specifically, for the omnibus statistic (Sk, Ku) 13-by-13 = (0.48, 2.95) and (Sk, Ku) 17-by-17 = (0.25, 2.98). These values are typically worse for the omnibus version of the statistic since the dependence structure in the multi-image data is more complex.

VIII. DISCUSSION AND CONCLUSIONS

In this paper we have proposed a novel statistical approach to patch-based change detection on two or more coregistered images. This approach gives a unified statistical thresholding procedure to perform change detection based on statistical features that have a known distribution under the no-change hypothesis. The proposed FDRformulation is based on the control of the proportion of false alarms among all detections. This efficient technique for large scale hypotheses testing allows to use the wide range of statistical tests developed in the state-of-the-art [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF] by adjusting to the dependence structure present in the images and the patch-based samples. The approach involves only a few parameters and is highly parallelizable. We have presented several rank-based statistical features that report accurate experimental results and the corresponding detectors positively compare with benchmark techniques in three different applications. Further features can be easily constructed to elaborate application-specific change detectors.

The developed change detection approach constitutes a flexible and statistically sound inference procedure that, nevertheless, has several limitations. First, we assume that the amount of changes present on images is limited (at most 25% of the total image size). This restriction is a price to pay for performing the change detection in an unsupervised manner with no learning data available (such as earlier frames). For the same reason the developed method underperforms in case of more complicated background motion or alteration scenario, e.g., waves or foliage motion for still camera imagery, where the temporal approaches [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF], [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF] are expected to perform well. This is due to the fact that in an unsupervised scenario analysing two or several images it is very challenging to tune the algorithm to disregard such false positives. The same limitation applies also to methods like [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF]. Second, the considered statistical features should have a known distribution under the no-change hypothesis. The presented experiments demonstrate that these limitations are not critical for various applications.

We have investigated the use of FDR-based change detection in case of multiple input images. To this end we have employed an omnibus Levene statistic to test the equality of variances. This statistic can also be employed in the two-sample scenario, and can be easily modified to test the differences in the means by suppressing step 2 in Algorithm 3. The omnibus test can be followed by post-hoc pairwise Levene tests to identify the specific image pairs that report meaningful changes. It is important to mention that post-hoc testing corresponds to a two-step procedure that may not ensure a strict control over the type I error. Furthermore, it is possible to have a change reported by omnibus statistic and no significant change witnessed by any pairwise test. The opposite could also occur.

The reason for that is the difference in the forms of the rejection regions specified by these two test procedures.

Nevertheless, a strong advantage of the omnibus statistic is the severely reduced computational complexity, and for some applications the specification of the time of the change is not relevant.

In all of the performed experiments we have observed only minor deviations of z-scores from the normal distribution as a result of intra-sample dependence, see Fig. 1 as an example, which experimentally validates the assumptions of Section IV. The good results obtained with different types of images in diverse applications are interpreted as due to the generality of the proposed approach. Note that the developed method (as well as similar spatial approaches [START_REF] Singh | Order consistent change detection via fast statistical significance testing[END_REF], [START_REF] Lanza | Statistical change detection by the pool adjacent violators algorithm[END_REF], [START_REF] Xie | Sudden illumination change detection using order consistency[END_REF]) underperforms in case of complicated background motion, when the temporal approaches [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF], [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF] perform better in discriminating back-and foreground of the images.

The main directions for the future research include: i) the integration of the developed change detector into a multi-resolution approach to automatically select the patch sizes, and (ii) incorporation of the patches' dependence structure into the model, for instance, via hidden Markov models.

The quantities ν jk = N i=1 (Y ij -Ȳ•j )(Y ik -Ȳ•k )/(N -1) for j = 1, . . . , K and k = 1, . . . , K are the clipped sample covariances. The following averages are employed above: overall ν = 1/K 2 j k ν jk , diagonal νd = 1/K j ν jj , and row-wise νj• = 1/K k ν jk .

The Levene statistic's z-values are calculated as z = Φ -1 (p), with the corresponding p-values computed through the Fisher distribution as p = 1 -F v1,v2 (f). The complete process is presented in Algorithm 3.

Fig. 1 .

 1 Fig. 1. The normalized histogram and its estimates for the Wilcoxon feature on the XSAR image pair. The z-scores in boxes correspond to the detection areas with FDR 0.1.

(

  temporal dependence or inter-sample dependence). Therefore, when considering the two-sample hypothesis tests on local patch samples we have to deal with two kinds of dependences in samples. For the problem of image change detection we propose to use the patch-based formulation of the FDR-based technique presented in Section III. The resulting empirical-Bayesian change detector allows to address: (i) the inference based on correlated individual tests (due to overlap of patch-based samples), and (ii) intra-sample dependence of patches.

0 :D

 0 X ≡ D Y }, where D X and D Y are the CDFs of the random variables X and Y observed by samples X and Y of size N , respectively. The CvM test can detect any kind of distribution deviations given sufficient sample sizes.

Algorithm 2 :

 2 Pairwise post-hoc testing 1 perform C = K(K -1)/2 pairwise Levene tests to obtain the corresponding p-values p 1 , . . . , p C ; 2 order the results so that p (1) . . . p (C) ; 3 set k = 0; 4 increment k by 1; 5 if k > C, stop; 6 if p (k) γ(C -k + 1)/C, reject hypotheses reporting p (k) , . . . , p (C) and stop, else goto step 4.(a) XSAR image 1 (b) XSAR image 2 (c) Ln-GKIT [25] (d) thresholded KL [31], 9 × 9 (e) FDR-based CvM, 9 × 9 (f) FDR-based Wilcoxon, 13 × 13

Fig. 2 .

 2 Fig. 2. (a),(b) 700 × 300 pixel XSAR image pair of Pavia, Italy, the largest changed areas are highlighted by red circles in (b). Various detection results (white -no change, black -change) in (c)-(f).

Fig. 3 .

 3 Fig. 3. ROCs for the employed detectors on XSAR and still camera datasets.

7 Fig. 4 . 5 (Fig. 5 .

 7455 Fig. 4. 650 × 200 pixel COSMO-SkyMed 2.5m resolution StripMap image pair ( c ASI) of Port-au-Prince, Haiti [(a),(b)], and various detection results (white -no change, black -change) in (c)-(h).

Fig. 6 .

 6 Fig. 6. (a),(b),(e) 360 × 240 and (c),(d) 720 × 580 pixels still camera image pairs from the database [16], and the change detection results by methods [9], [24] and the FDR-based approach with (a)-(c) mCvM, (d) Wilcoxon and (e) CvM features.

Fig. 7 .

 7 Fig. 7. Semi-synthetic SAR image sequence: (a) ground truth and (b)-(d) images, results of (e),(f) pairwise and (g),(h) omnibus Levene FDR-detection.

TABLE I XSAR

 I PAIR DETECTION RESULTS: FALSE POSITIVE RATES (FPR), TRUE POSITIVE RATES (TPR), FALSE DISCOVERY PROPORTIONS (FDP),

	TOTAL AMOUNTS OF DETECTION AND COMPUTATION TIMES (SECONDS)
	Method	FPR	TPR	FDP	Detection	Time
	Ln-GKIT [25]	0.58% 90.36% 9.78%	3.24%	1.08
	KL-based [31]	0.40% 86.67% 7.12%	3.10%	2.03
	FDR CvM	0.08% 94.76% 1.33%	3.01%	1.54
	FDR Wilcoxon 0.28% 98.80% 4.55%	3.69%	2.84
	strategy could consist in pooling all the pairwise tests from all pixels rejected by the omnibus Levene FDR-analysis,
	and performing a second round of FDR-detection targeting image pairs according to Algorithm 2.

ACKNOWLEDGMENT

The authors would like to thank Prof. P. Gamba (University of Pavia, Italy) for providing the XSAR dataset, and the Italian Space Agency for providing the COSMO-SkyMed (CSK R ) dataset of Haiti (COSMO-SkyMed Product c ASI -Agenzia Spaziale Italiana -2010. All Rights Reserved).

This work was partially supported by the Italian Space Agency (ASI) in the framework of the project entitled "Development and

next apply the clipping that allows to improve robustness of the statistic to ouliers. Specifically, let

where O R kn is the ordered version of O kn such that values of O R kn in each group k are in ascending order, [ωN ] is the integer part of ωN , and ω is the amount of clipping. The modulus transform employed in O kn projects the "extreme" observations on the right side, and hence the clipping is also one-sided.

The Levene statistic for the dependent samples [START_REF] Efron | Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction[END_REF] requires

Here ŌT k denotes the mean of the ω-trimmed sample k (i.e. calculated by ignoring the clipped observations), Ō = ŌT k /K, and averages

The critical values for the statistic F are taken from the Fisher CDF F v1,v2 with degrees of freedom v 1 = ε(K -1)

. Compared to the textbook ANOVA [START_REF] Conover | Practical nonparametric statistics[END_REF], the ε-correction allows to relax the technical assumptions of sphericity of the variance and correlation structure of the K samples, see [START_REF] Wilcox | Applying contemporary statistical techniques[END_REF]. The quantity ε is estimated as ε = N (K -1)A/B -2 (K -1)(N -1 -(K -1)A/B) ,


Christophe Alias

Carsten Fuhs

Laure Gonnord

Estimation of Parallel Complexity with Rewriting Techniques

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The motivation of this work is the automatic transformation of sequential code into parallel code without changing its (big-step operational) semantics, only changing the computation order is allowed. We want to find out the limits of these approaches, by characterising the "maximum level of parallelism" that we can find for a given sequential implementation of an algorithm.

In this paper, we propose a way to estimate the parallel complexity which can be informally defined as the complexity of the program if it were executed on an machine with unbounded parallelism.

Such a result could come as by-product of the polyhedral-based automatic parallelisation techniques for array-based programs [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, I, one-dimensional time[END_REF]. However, for general programs with complex data flow and inductive structures, such techniques have not been explored so far.

The contribution of this paper is a novel unifying way to express program dependencies for general programs with inductive data structures (lists, trees. . . ) as well as a way to use complexity bounds of term rewriting systems in order to derive an estimation of this parallel complexity.

2 Intuitions behind the notion of "parallel complexity" for ( i =0; i <= N ; i ++) for ( j =0; j <= N ; j ++)

// Block S { m1 [ i ][ j ] = Integer . MIN_VALUE ; for ( k =1; k <= i ; k ++) m1 [ i ][ j ] = max ( m1 [ i ][ j ] , H [i -k ][ j ] + W [ k ]); m2 [ i ][ j ] = Integer . MIN_VALUE ; for ( k =1; k <= j ; k ++) m2 [ i ][ j ] = max ( m2 [ i ][ j ] , H [ i ][ j -k ] + W [ k ]); H [ i ][ j ] = max (0 , H (i -1 ,j -1)+ s ( a [ i ] , b [ i ]) , m1 [ i ][ j ] , m2 [ i ][ j ]); } j i Figure 1:
The Smith-Waterman sequence alignment algorithm and its dependencies. Each point (i, j) represents an execution of the block S, denoted by S, i, j . In this paper, we consider the derivation of a lower bound for the intrinsic parallelism of a sequential imperative program P and thus an upper bound for the complexity of a fully parallelised version of P . Figures 1 and2 depict the two motivating examples that will be studied in the remainder of this paper. Figure 1 is a loop kernel computing the Smith-Waterman optimal sequence alignment algorithm1 . Figure 2 is a simple recursive function to compute the maximum element of a binary tree. Usually, the parallelism is found by analysing the data dependencies between the operations executed by the program. There is a data dependency from operation o 1 to operation o 2 in the execution of a program if o 1 is executed before o 2 and both operations access the same memory location. In particular, we have a flow dependency when the result of o 1 is required by o 2 . Non-flow dependencies (write after write, write after read) can be removed by playing on the memory allocation [START_REF] Barthou | Maximal static expansion[END_REF] and can thus be ignored. Thus, we consider only flow dependencies referred as dependencies in the remainder of this paper. If there is no dependency between two operations, then they may be executed in parallel. While the dependency relation is in general undecidable, in practice we can use decidable over-approximations such that a statement that two operations are independent is always correct.

In Figure 1, each point represents an execution of the block S computing H[i][j] for a given i and j in 0, N . Such an operation is written S, i, j . The arrows represents the dependencies towards a given S, i, j . For instance the diagonal arrow means that

H[i][j] requires the value of H[i-1][j-1] to be computed.
In Figure 2, we depicted the execution trace of the recursive program on a tree. Here, the dependencies between computations resemble the recursive calls: the dependency graph (in the sense used for imperative programs) is the call tree.

All reorderings of computations respecting the dependencies are valid orderings (that we will name schedule in the following). In both cases, the dependencies we draw show a certain potential parallelism. Indeed, each pair of computations that do not transitively depend on each other can be executed in parallel; moreover, even a machine with infinite memory and infinitely parallel cannot do the computation in an amount of time which is shorter than the longest path in the dependency graph. The length of this longest path, referred to as parallel complexity, is thus an estimation of the potential parallelism of the program (its execution time with suitably reordered instructions on an idealised parallel machine).

Our goal in this paper is, given a sequential program P and an over-approximation of its dependency relation, to find bounds on the parallel complexity of P and the over-approximation of its dependency relation. Via the representation of the dependencies via (possibly constrained) rewrite rules, we can apply existing techniques to find such bounds for (constrained) rewrite systems (e.g., [START_REF] Hirokawa | Automated complexity analysis based on the dependency pair method[END_REF][START_REF] Brockschmidt | Analyzing runtime and size complexity of integer programs[END_REF]).

Computing the parallel complexity of programs

In this section we explain on the two running examples the relationship between termination and scheduling (this relationship was already explored a bit in [START_REF] Alias | Multi-dimensional rankings, program termination, and complexity bounds of flowchart programs[END_REF]), and polynomial interpretations (more generally, proofs of complexity bounds for term rewriting) and parallel complexity. As the long-term goal is to devise compiler optimisations by automatic parallelisation for imperative programs, we consider programs with data structures such as arrays, structs, records and even classes. The main idea is that all these constructions can be classically represented as terms via the notion of position, and the dependencies as term rewriting rules. Contrarily to other approaches for proving termination of programs by an abstraction to rewriting (e.g., [START_REF] Otto | Automated termination analysis of Java Bytecode by term rewriting[END_REF][START_REF] Falke | Termination analysis of C programs using compiler intermediate languages[END_REF]), we only encode the dependencies (and forget about the control flow).

First example: Smith-Waterman algorithm

For the Smith-Waterman program of Figure 1, we can derive (with polyhedral array dataflow analysis, as in [START_REF] Barthou | Maximal static expansion[END_REF]), forgetting about the local computations of W scores, the following dependencies as constrained rewrite rules:

dep(i, j) → dep(i -1, j -1) : 0 ≤ i ≤ n, 0 ≤ j ≤ n dep(i, j) → dep(i -k, j) : 0 ≤ i ≤ n, 0 ≤ j ≤ n, 1 ≤ k ≤ i dep(i, j) → dep(i, j -) : 0 ≤ i ≤ n, 0 ≤ j ≤ n, 1 ≤ ≤ j
To get an estimation of the parallel complexity, following the inspiration of previous work on termination proofs for complexity analysis (e.g., [START_REF] Hirokawa | Automated complexity analysis based on the dependency pair method[END_REF][START_REF] Alias | Multi-dimensional rankings, program termination, and complexity bounds of flowchart programs[END_REF][START_REF] Brockschmidt | Analyzing runtime and size complexity of integer programs[END_REF]), we first try to generate a polynomial interpretation [START_REF] Lankford | On proving term rewriting systems are Noetherian[END_REF] to map function symbols and terms to polynomials:

Pol (dep(x 1 , x 2 )) = x 1 + x 2
Essentially, this says that iteration S, x 1 , x 2 can be executed at time stamp x 1 + x 2 . As 0 ≤ x 1 , x 2 ≤ N , the maximal timestamp is 2N = O(N ). Not only does this mean that the program may be parallelised, but it provides an actual reordering of the computation, along the parallel wavefront

x 1 + x 2 .
This linear interpretation (or ranking function) provides us with a bound in O(N 1 ) for the parallel complexity of the program. Thus, the degree of sequentiality is 1. As the overall runtime of the program is in O(N 2 ), this gives us a degree of parallelism of 2 -1 = 1 [START_REF] Karp | The organization of computations for uniform recurrence equations[END_REF].

Let us point out that we would have obtained the same results using affine scheduling techniques from the polyhedral model [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, I, one-dimensional time[END_REF]. The interesting fact here is that our apparatus is not restricted to regular programs (for loops, arrays) as the polyhedral model. Also, current complexity analysis tools like KoAT [START_REF] Brockschmidt | Analyzing runtime and size complexity of integer programs[END_REF] are able to compute similar results within a reasonable amount of time. The next section show how our technique applies to a recursive program on inductive data structures. An object of class Tree is represented by the term Tree(val , left, right), for some terms val , left, right that represent its attributes [START_REF] Otto | Automated termination analysis of Java Bytecode by term rewriting[END_REF]. For instance, the Java object defined in Figure 3 corresponds to the following term:

Second example: computing the maximum element in a tree

t = Tree( 2, Tree(3, null, null), Tree(4, Tree(7, null, null), null))
Recall that to address entries of an n-dimensional array, we use vectors (i 1 , . . . , i n ) ∈ N n as indices or "positions". Then for an array A, we say that A[i 1 ]. . .[i n ] is the "entry" at the array's position (i 1 , . . . , i n ). Similarly to array entries, we would also like to address particular "entries" of a term, i.e., its subterms.

For more general terms, we can use a similar notion: positions. For our tree t, we have Pos(t) = {ε, 1, 2, 3, 21, 22, 23, . . .}, giving an absolute way to access each element or subarray. For instance, 21 denotes the first element of the left child of the tree (i.e., the number 3). Now let us consider the program in Figure 2. This function computes recursively the maximum value of an integer binary tree. Clearly, the computation of the maximum of a tree depends on the computation of each of its children. However, the computation of the max of each child is independent from the other. There is thus potential parallelism.

Here we observe structural dependencies from the accesses to the children of the current node. Like in the previous example, from the program we generate the following "dependencyrewrite rules" (note that similar to the dependency pair setting for termination proving [START_REF] Arts | Termination of term rewriting using dependency pairs[END_REF], it suffices to consider "top-rewriting" with rewrite steps only at the root of the term): We can use the following polynomial interpretation (analogous to the notion of a ranking function) to prove termination and also a complexity bound:

Pol (dep(x 1 )) = x 1 and Pol (Tree(x 1 , x 2 , x 3 )) = x 2 + x 3 + 1
Or, using interpretations also involving the maximum function [START_REF] Fuhs | Maximal termination[END_REF]:

Pol (dep(x 1 )) = x 1 and Pol (Tree(x 1 , x 2 , x 3 )) = max(x 2 , x 3 ) + 1
Thus we interpret a tree as the maximum of its two children plus one to prove termination of the dependency relation of the original sequential program, essentially mapping a tree to its height. This means that the parallel complexity (i.e., the length of a chain in the dependency relation) of the program is bounded by the height of the input data structure. Indeed, the two recursive calls could be executed in parallel, with a runtime bounded by the height of the tree on a machine with unbounded parallelism. The interpretation Pol (Tree(x 1 , x 2 , x 3 )) = max(x 2 , x 3 ) + 1 induces a wavefront for the parallel execution along the levels of the tree, i.e., the nodes at the same depth in the tree.

In contrast, the overall runtime of the original sequential program is bounded only by the size of the input tree, which may be exponentially larger than its depth.

Conclusion and Future Work

In this paper we showed some preliminary results on the (automatable) computation of the parallel complexity of programs with inductive data structures.

In the future, we will investigate a complete formalisation of these preliminary results, and test for applicability in more challenging programs like heapsort and prefixsum. As we said in the introduction, expressing the parallel complexity is the first step toward more ambitious use of rewriting techniques for program optimisation. The work in progress includes the computation of parallel schedules from the rewriting rules or their interpretation, and then parallel code generation from the obtained schedules.

Figure 2 :

 2 Figure 2: Maximum element of a tree algorithm and its call tree.

Figure 3 :

 3 Figure 3: The tree t.

  dep(Tree(val , left, right)) → dep(left) (S1) dep(Tree(val , left, right)) → dep(right) (S2)

See https://en.wikipedia.org/wiki/Smith-Waterman_algorithm. We consider two sequences of the same length N .

* This work was partially supported by a "BQR" funding at ENS De Lyon.