Abdul-Fattah Abu-Awwad

Véronique Maume-Deschamps

Pierre Ribereau

Abu-Awwad Abdul-Fattah

Maume-Deschamps V Éronique

Ribereau Pierre

Fitting spatial max

Keywords: Extremal dependence, Madogram, Max-stable model, Max-mixture model, Nonlinear least squares

published or not. The documents may come

Introduction

Max-stable stochastic processes arise as a fundamental class of models that are able to describe spatial extreme value phenomena. Max-stable process models for spatial data were first constructed using the spectral representation [START_REF] De | A spectral representation for max-stable processes[END_REF]. Several subsequent works on the construction of spatial max-stable processes have been developed, see e.g. [START_REF] Richard | Max-stable processes and spatial extremes[END_REF][START_REF] Schlather | Models for stationary max-stable random fields[END_REF][START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF][START_REF] Anthony | [END_REF]. The inference on spatial processes is an open field that is still in development. Both parametric and nonparametric inference methods are used in the literature.

Camille Jordan Institute -ICJ, Université Claude Bernard Lyon 1, France Address for correspondence: A. Abu-Awwad, Department of Applied Mathematics, School of informatics and mathematics, Camille Jordan Institute -ICJ, Université Claude Bernard Lyon 1, 43 boulevard du 11 novembre 1918, F-69622 Villeurbanne Cedex, France. E-mail addresses: abuawwad@math.univ-lyon1.fr, veronique.maume@univ-lyon1.fr, pierre.ribereau@univ-lyon1.fr.

For a stationary spatial process X = {X(s), s ∈ S}, the Asymptotic Dependence (AD) is characterized by χ(h) > 0, with:

(1.1)

χ(h) = lim x→∞ P{X(s) > x|X(s + h) > x}, s, s + h ∈ S.
This means that, for an AD process, a large event at location s + h leads to a non-zero probability of a similarly large event at location s for some spatial lag vector h. On one other hand, a process is Asymptotically Independent (AI) if χ(h) = 0 for any h. This is achieved e.g. for Gaussian processes, see [START_REF] Sibuya | Bivariate extreme statistics, i[END_REF].

Within the class of max-stable models, only two types of dependence structures are possible: either the process is AD or it is independent. This restriction leads to a drawback of max-stable processes that they are too coarse to describe multivariate joint tails with asymptotic independence sufficiently accurately. Particularly, fitting asymptotic dependent models to asymptotically independent data leads over/under estimation of probabilities of extreme joint events, since there is a mis-placed assumption that the most extreme marginal events may occur simultaneously [START_REF] Coles | Dependence measures for extreme value analyses[END_REF]. [START_REF] Thibaud | Threshold modeling of extreme spatial rainfall[END_REF]Davison et al., 2013) introduced recent examples about practical difficulties to identify whether a data set should be modeled using an asymptotically dependent or asymptotically independent. [START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF] introduced a new class of models, so-called max-mixture models, to capture both AD and AI. The basic idea is to mix max-stable and asymptotic independent processes. Let a ∈ [0, 1], then the max-mixture (MM) model is defined as

Z(s) = max{ aX(s), (1 -a)Y (s)}, s ∈ S
where X(s) is a stationary max-stable process and Y (s) is a stationary AI process, so that the parameter a represents the proportion of AD in the process Z.

In this paper, we are concerned with constructing a model-free criterion to choose a realistic value for the mixing parameter a. Our objective is not to model extremal dependence of joint tails but to set up a statistical criterion that facilitate the modelling of the spatial data with suitable behaviour. We shall use least squares on the F λ -madogram. In [START_REF] Bel | Assessing extremal dependence of environmental spatial fields[END_REF], a madogram based test on the AD of a max-stable process is proposed, while in [START_REF] Abu-Awwad | Censored pairwise likelihood-based tests for mixing coefficient of spatial max-mixture models[END_REF], a parametric test on a for max-mixture processes is developed.

The paper is organized as follows. Section 2 reviews spatial extremes processes. The proposed F λ -madogram for max-mixture models and the selection criterion for the mixing coefficient a are developed in Section 3, while Section 4 illustrates the performance of our method through a number of simulation studies. We conclude with an illustration of spatial analysis of precipitation in Australia in Section 5.

Spatial extremes processes: models

Throughout our work, X = {X(s), s ∈ S}, S ⊂ R d (generally, d = 2) is a spatial process, it will be assumed to be stationary and isotropic.

2.1. Max-stable processes. Suppose that {Y i (s) :

s ∈ S ⊂ R d }, i = 0, 1, 2, ..., are i.i.d
replicates of a random process Y (s), and that there are sequences of continuous functions {a n (s) > 0} and {b n (s)} such that, the rescaled process of maxima,

(2.1)

n i=1 Y i (s) -b n (s) a n (s) d → X(s), n → ∞
where the limiting random process X is assumed to be non-degenerate. By [START_REF] De | Spatial extremes: Models for the stationary case[END_REF] the class of the limiting processes X(s) coincides with the class of max-stable processes. This definition of MS processes offers a natural choice for modeling spatial extremes.

The univariate extreme value theory, implies that the marginal distributions of X(s) are Generalized Extreme value (GEV) distributed, and without loss of generality the margins can transformed to a simple MS process called standard Fréchet distribution, P(X(s) ≤ z) = exp{-z -1 }.

Following [START_REF] De | A spectral representation for max-stable processes[END_REF][START_REF] Schlather | Models for stationary max-stable random fields[END_REF], a simple MS process X(s) has the following representation (2.2)

X(s) = max k≥1 Q k (s)/P k , s ∈ S.
where Q k (s) are independent replicates of a non-negative stochastic process Q(s) with unit mean at each s, and P k are the points of a unit rate Poisson process (0, ∞).

For K ∈ N \ {0}, s 1 , . . . , s K ∈ S, and x 1 , . . . , x K > 0, the finite K-dimensional distribution of the process X owing to the representation( 2.2) is given by

(2.3) -log P(X(s 1 ) ≤ z 1 , . . . , X(s K ) ≤ z K ) = E K k=1 Q(s k ) z k = V s 1 ,...,s K (z 1 , . . . , z K )
where V s 1 ,...,s K (.) is called the exponent measure. It summarises the structure of extremal dependence, and satisfies the property of homogenity of order -1 and V s 1 ,...,s K (∞, ..., z, ..., ∞) = z -1 . It has to be noted that

-z log P{X(s 1 ) ≤ z, .., X(s K ) ≤ z} = V s 1 ,...,s K (1, ..., 1) = θ s 1 ,...,s K ,
The coefficient θ s 1 ,...,s K is known as the extremal coefficient. It can be seen as a summary of extremal dependence with two boundary values. The complete independence is achieved when θ s 1 ,...,s K = 1, while complete independence is achieved when θ s 1 ,...,s K = K. In the bivariate case, the AI and AD between a pair of random variables Z 1 and Z 2 , with marginal distributions F 1 and F 2 , may be identified by

(2.4) χ = lim u → 1 - P (F 1 (Z 1 ) > u|F 2 (Z 2 ) > u)).
The cases χ = 0 and χ > 0 represent AI and AD, respectively, [START_REF] Joe | Parametric families of multivariate distributions with given margins[END_REF]. This coefficient is related to the pairwise extremal coefficient θ through the relation χ = 2 -θ.

Since both dependence functions θ and χ are useless for AI processes, [START_REF] Coles | Dependence measures for extreme value analyses[END_REF] proposed a new dependence coefficient which measures the strength of dependence for AI processes:

(2.5) χ = lim

u → 1 - χ(u) = lim u → 1 - 2 log P (F (Z(s)) > u) log P (F (Z(s)) > u, F (Z(s + h)) > u) -1 AD (respectively AI) is achieved if and if χ = 1 (resp. χ < 1).
Another dependence model for bivariate joint tails was introduced by (Ledford and Tawn, 1996)

(2.6) P(X(s) > x, X(s + h) > x) ∼ x -1/η(h) L h (x), for x → ∞
where X is a stationary spatial process with unit Fréchet margins, L h (.) is a slowly varying function at ∞ and the tail dependence coefficient η(h) ∈ (0, 1]. AI corresponds to η(h) < 1.

Different choices for the process Q(s) in (2.2) lead to more or less flexible models for spatial maxima. Commonly used models are the Guassian extreme value process [START_REF] Richard | Max-stable processes and spatial extremes[END_REF], the extremal Gaussian process [START_REF] Schlather | Models for stationary max-stable random fields[END_REF], the Brown-Resnick process [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF], and the extremal-t process [START_REF] Opitz | Extremal t processes: Elliptical domain of attraction and a spectral representation[END_REF]. Below, we list these four specific examples of max-stable models.

The storm profile model [START_REF] Richard | Max-stable processes and spatial extremes[END_REF], is defined by taking

Q k (s) = f (s -W k ) in Equation (2.2)
, where f is the density function of a Gaussian random vector with covariance matrix Σ ∈ R 2×2 . The function f plays a major role as it determines the shape of the storm events.

W k is a homogenous Poisson process. The bivariate exponent function of the Smith model has the form

-log P[X(s) ≤ x 1 , X(s + h) ≤ x 2 ] = 1 x 1 Φ β(h) 2 + 1 β(h) log x 2 x 1 + 1 x 2 Φ β(h) 2 + 1 β(h) log x 1 x 2
where β(h) = √ h T Σ -1 h and Φ is the standard normal distribution function. In this case the extremal coefficient is equal to θ(h) = 2Φ{β(h)/2}.

A model originally due to [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] is the Truncated Extremal Gaussian (TEG) model and has been exemplified in [START_REF] Anthony | [END_REF]. This process is ob- The bivariate exponent function of a TEG model in the stationary case has the form (2.7)

tained by taking Q k (s) = c max(0, ε k (s))1 A k (s -R k ),
-log P[X(s) ≤ x 1 , X(s + h) ≤ x 2 ] = 1 x 1 + 1 x 2 1 - α(h) 2 1 -1 - 2(ρ(h) + 1)x 1 x 2 (x 1 + x 2 ) 2 where α(h) = (1 -h/2r)1 [0,2r] if A is a disk of fixed radius r. This yields θ(h) = 2 -α(h){1 - (1 -ρ(h))/2}.
The max-stable Brown-Resnick (BR) process model proposed by [START_REF] Bruce | Extreme values of independent stochastic processes[END_REF][START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]) is a stationary max-stable process that can be constructed with Q k (s) = exp{ε k (s) -γ(s)}, s ∈ S, where ε k (s) denotes a Gaussian process with semivariogram γ(h). The bivariate exponent function of a BR process is:

-log P[X(s) ≤ x 1 , X(s + h) ≤ x 2 ] = 1 x 1 Φ γ(h) 2 + 1 2γ(h) log x 2 x 1 + 1 x 2 Φ γ(h) 2 + 1 2γ(h) log x 1 x 2
where γ and Φ denote respectively the semivariogram and the standard normal distribution function. In particular, when the variogram 2γ(h) = h T Σ -1 h for some covariance matrix Σ, we recover the bivariate distribution function of a Smith model.The pairwise extremal coefficient

for a Brown-Resnick process is θ(h) = 2Φ{ γ(h)/2}.
The extremal-t max-stable process proposed in [START_REF] Opitz | Extremal t processes: Elliptical domain of attraction and a spectral representation[END_REF][START_REF] Ribatet | Extreme value copulas and max-stable processes[END_REF]) can be constructed by using

Q k (s) = m -1/v v T k (s) v
, where T k is a zero mean Gaussian process

with correlation function ρ, v ≥ 1, 1/m v = √ π2 v/2-1 Γ( v+1 2
), Γ(.) is the gamma function. This process has the bivariate exponent function:

-log P[X(s) ≤ x 1 , X(s + h) ≤ x 2 ] = 1 x 1 T v+1 αρ(h) + α x 2 x 1 1/v + 1 x 2 T v+1 αρ(h) + α x 1 x 2 1/v
where T v is the distribution function of a Student random variable with v degrees of freedom

and α = [v + 1/{1 -ρ 2 (h)}] 1/2
. For an extremal-t process the degrees of freedom v controls the upper bound of the extremal coefficient:

θ(h) = 2T v+1 (v + 1) 1-ρ(h) 1+ρ(h) .
In this paper, we shall make intensive use of the so-called F -madogram that is based on a classical geostatistical tool, the madogram [START_REF] Matheron | Suffit-il, pour une covariance, dâêtre de type positif[END_REF] . It has been introduced in [START_REF] Cooley | Variograms for spatial max-stable random fields[END_REF]. Let X is a stationary max-stable random process. The marginal distribution function is denoted by F . The F -madogram is defined by:

(2.8)

ν F (h) = 1 2 E[|F (X(s)) -F (X(s + h))|], 0 ≤ ν F (h) ≤ 1/6.
The bounds of the F madogram correspond respectively to complete dependence and independence. Due to the one-to-one relationship between the extremal dependence function and the F -madogram, a simple estimator for θ(h) can be derived:

(2.9)

θ(h) = 0.5 + ν F (h) 0.5 -ν F (h)
where

ν(h) = 1 2N N i=1 | F {x i (s)} -F {x i (s + h)}|, x i (
s) and x i (s + h) are the i-th observations of the random field at locations s and s + h. F is the empirical distribution function, i.e,

F (z) = (N + 1) -1 N i=1 I Z i (s j )≤z
, where I(.) is the indicator function.

The so-called F λ -madogram has been introduced for max-stable models by [START_REF] Bel | Assessing extremal dependence of environmental spatial fields[END_REF] as a generalization of F -madogram (2.8): for any λ > 0, let

(2.10) ν F λ (h) = 1 2 E |F λ {X(s)} -F λ {X(s + h)}| .
A nonlinear least squares procedure has been proposed by [START_REF] Bel | Assessing extremal dependence of environmental spatial fields[END_REF] based on F λ -madogram to compute an estimator for the extremal dependence function that may outperforms other known estimators. In that work, it has been found by some trials that good estimations are obtained for λ ∈ [2, 3].

2.2. Hybrid models of spatial extremal dependence. Although max-stable models seem to be suitable for modeling extremely high threshold exceedances, asymptotic independence models may show a better fit at finite thresholds. Due to difficulty or impossibility in practice to decide whether a dataset should be modeled using AD or AI, [START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF] have been introduced the hybrid spatial dependence models which are able to capture both AD and AI.

Consider Y a stationary Gaussian process. Let Y (s) = -1/ log(Φ(Y (s))) then, Y is an AI process with unit Frchet marginal distributions. Another class of AI processes called inverted max-stable processes has been proposed by [START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF]. They are defined as

(2.11) Y (s) = -1/ log{1 -exp [-Y (s) -1 ]}
where Y is a simple max-stable process with extremal coefficient θ Y . We a slight abuse of language, we shall denote θ Y by θ Y . With this construction, each max-stable process may be transformed into an AI independent counterpart. This inverted max-stable process (IMS) satisfies (2.6) and η(h) = 1/θ Y (h). The bivariate distribution function is given by

P(Y (s) ≤ y 1 , Y (s + h) ≤ y 2 ) = -1 + exp(-y -1 1 ) + exp(-y -1 2 ) + exp{-V Y [ω(y 1 ), ω(y 2 )]} where V Y is the exponent measure of the bivariate extreme-distribution of {Y (s), Y (s + h)},
and

ω(y) = -1/ log{1 -exp [-y -1 ]}.
We a slight abuse of language, we shall say that V Y is the exponent measure of Y .

We are now in position to define the max-mixture processes that we will be working on. Let X be a simple max-stable process with bivariate extremal coefficient θ X , and Y be an inverted max-stable process with coefficient of tail dependence η. Assume that X and Y are independent. Then for a mixture proportion a ∈ [0, 1], the spatial max-mixture process proposed by [START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF] is defined as

(2.12) Z(s) = max{ aX(s), (1 -a)Y (s)}.
Clearly, models that are only AD or AI are submodels of Z, obtained for a = 1, a = 0, respectively. The bivariate joint survivor function of the process Z satisfies

(2.13) P[Z(s) > z, Z(s + h) > z] = a(2 -θ(h)) z + 1 -a z 1/η(h) + O(z -2 ) as z → ∞ If h 0 = inf {h : θ(h) = 2}
is finite, then the process Z is AD up to distance h 0 , and AI for longer distances. The bivariate CDF for a pair of sites (Z(s), Z(s + h)) is straightforwardly obtained by the independence between X(s), Y (s)

(2.14)

P(Z(s) ≤ z 1 , Z(s+h) ≤ z 2 ) = P X(s) ≤ z 1 a , X(s + h) ≤ z 2 a P Y (s) ≤ z 1 1 -a , Y (s + h) ≤ z 2 1 -a
Thus, in the case where X(s) is a max-stable process and Y (s) is a inverted max-stable process, the distribution function in (2.14) has the form

P(Z(s) ≤ z 1 , Z(s + h) ≤ z 2 ) = exp{-aV X (z 1 , z 2 )} × {-1 + exp[(a -1)/z 1 ] (2.15) + exp[(a -1)/z 2 ] + exp[-V Y [ω((1 -a)/z 1 ), ω((1 -a)/z 2 )]}
where V X and V Y are the bivariate exponent measures for X and Y respectively.

Figure 1 displays two simulated images of AI processes over the [0, 10] 2 square. The corresponding functions η(h) are also represented. According to [START_REF] Anthony | Modelling dependence within joint tail regions[END_REF], the case η(h) = 1/2 corresponds to the near-independence, the AI process constructed from a Brown-Resnick process (this case corresponds to the isotropic Smith process) allows asymptotic independence but tends to near-indpendence for long distances. While an AI process constructed from extremal-t process presents a stronger dependence in the asymptotic independence when h is sufficiently large. (2.13), the measure χ Z associated to the process Z can be computed for a distance h as χ Z (h) = aχ X (h), see [START_REF] Bacro | A flexible dependence model for spatial extremes[END_REF]. The hybrid model extends traditional dependence modeling within the AD class and is appropriate when AD is present at all distances because it permits to capture a second order in the dependence structure which is not possible with a max-stable model. ) functions for a max-mixture process Z(.) based on two max-stable processes X(.) and X (.) that belong to the same family with the same parameters used in Figure 1 for a ∈ {0.25, 0.5, 0.75}.

In [START_REF] Bacro | A flexible dependence model for spatial extremes[END_REF] daily rainfall data in the East of Australia are studied. Different models (MS, AI, and MM) are fitted to the data. It is showed that MM models has the merit to overcome the limits of MS models in which only AD or exact AI can be modeled.

F λ -madogram for spatial max-mixture model

In the present paper we shall use the F λ -madogram defined by Equation (2.10) for maxmixture models (2.12). We begin by calculating the expression of the F λ -madogram for maxmixture models. Then, we shall develop a F λ -madogram procedure to estimate θ X , θ Y and chose a.

Proposition 3.1. Let {X(s), s ∈ S} be a simple max-stable process, with extremal coefficient function θ X (h), and Y (s) be an inverted max-stable process with coefficient of tail dependence

function η(h) = 1/θ Y (h). Let a ∈ [0, 1] and Z = max{aX, (1 -a)Y }. Then, the F λ -madogram
of the spatial max-mixture process Z(s) is given by

ν F λ (h) = λ 1 + λ - 2λ a(θ X (h) -1) + 1 + λ + λ aθ X (h) + λ - λθ Y (h) (1 -a)θ Y (h) + aθ X (h) + λ β aθ X (h) + λ 1 -a , θ Y (h) (3.1)
where β(., .) is the beta function.

As a consequence of Proposition 3.1, we easily recover the expressions of F λ -madograms for max-stable and inverted max-stable processes. The F λ -madogram of a simple max stable process X with extremal dependence coefficient θ X is:

(3.2) ν F λ (h) = λ λ + 1 θ X (h) -1 λ + θ X (h) and we have ν F λ (h) ∈ [0, λ (1+λ)(2+λ) ]. The F λ -madogram of an inverted max-stable process Y with extremal dependence coefficient θ Y is: (3.3) ν F λ (h) = 1 1 + λ - λθ Y (h) λ + θ Y (h) β(λ, θ Y (h))
Proof. of Proposition 3.1 We use that for any x, y ∈ R, |x -y| = 2 max {x, y} -(x + y) with

x = F λ {Z(s)} and y = F λ {Z(s + h)}. Moreover, recall that E[F α {Z(s)}] = 1/(1 + α). We have ν F λ (h) = E max{F λ {Z(s)}, F λ {Z(s + h)}} - 1 2 E F λ {Z(s)} + F λ {Z(s + h)} = E max{F λ {Z(s)}, F λ {Z(s + h)}} - 1 (1 + λ) Consider the random variable W = max{F λ {Z(s)}, F λ {Z(s + h)}}. Then probability distri- bution function G of W satisfies G(z) = P[W ≤ z] = P max{F λ {Z(s)}, F λ {Z(s + h)}} ≤ z = P Z(s) ≤ F -λ (z), Z(s + h) ≤ F -λ (z) = P Z(s) ≤ - λ log(z) , Z(s + h) ≤ - λ log(z) = exp -V X - λ a log(z)
, -λ a log(z)

• z

1-a λ + z 1-a λ -1 + exp    -V Y   - 1 log 1 -z 1-a λ , - 1 log 1 -z 1-a λ       
This rewrites:

G(z) = z aθ X (h) λ • 2z (1-a) λ -1 + 1 -z (1-a) λ θ Y (h) = 2z a(θ X (h)-1)+1 λ -z a λ θ X (h) + z a λ θ X (h) 1 -z (1-a) λ θ Y (h)
.

Thus, we are led to

E[W ] = 1 0 zdG = zG 1 0 - 1 0 Gdz = 1 - 2λ a(θ X (h) -1) + 1 + λ - λ aθ X (h) + λ + I and
I := 1 0 z a λ θ X (h) 1 -z (1-a) λ θ Y (h) dz = λ (1 -a) β aθ X (h) + λ 1 -a , θ Y (h) + 1 = λθ Y (h) (1 -a)θ Y (h) + aθ X (h) + λ β aθ X (h) + λ 1 -a , θ Y (h)
This proves Equation (3.1). The asymptotic independent process Y is an inverse TEG process with A Y a disk of fixed radius r Y = 0.5 and ρ Y (h) = exp(-h/0.4). Middle panel: Max-mixture model in which X is a TEG process as before.

The asymptotic independent process Y is an inverse extremal-t process with v = 1 degrees of freedom and ρ Y (h) = exp(-h/0.6). Right panel: Max-mixture model in which X is a TEG process as before.

Y is an inverted Brown-Resnick process with semivariogram γ(h) = h 2 . The mixing coefficient is setted to a = 0.5 in the three models. The grey vertical lines represent the diameters of the disks for TEG processes.

We are now in position to describe a choice scheme for the mixing parameter a of a maxmixture process Z in (2.12). From Equation (3.1), we may write the F λ -madogram as a function of a, λ, θ X and θ Y , that is ν F λ (h) = Φ(a, λ, θ X (h), θ Y (h)). The idea of our choice procedure is that θ X and θ Y may be estimated by θ X and θ Y , minimizing the square difference between Φ(a, λ, θ X (h), θ Y (h)) and its empirical counterpart, then we can choose a such that the empirical version of the F λ -madogram is the closest to Φ(a, λ , θ X (h), θ Y (h)). This idea is close to the non parametric estimation of the parameters of MM processes that has been proposed in [START_REF] Ahmed | A semiparametric estimation for max-mixture spatial processes[END_REF] as an alternative to maximum composite likelihood estimation.

Formally, we consider Z i , i = 1, ..., N copies of Z,

Q i (h, λ) = 1 2 |F λ (Z i (s)) -F λ (Z i (s + h))|,
where F denotes the unit Frchet distribution function. From the definition of the F λ -madogram,

we have E[Q i (h, λ)] = ν F λ (h).
Denote by Λ ⊂ [0, ∞) a finite set of some possible λ choices, then for a given value of a, a semi-parametric nonlinear least squares minimization procedure for estimating the extremal coefficient θ

(h) = (θ X (h), θ Y (h)) t is (3.4) θ a N LS (h) = arg min θ∈[1,2] 2 N -1 λ∈Λ i=1,...,N [Q i (h, λ) -Φ(a, λ, θ X (h), θ Y (h))] 2
Assume that the Z i 's are observed at locations s 1 , . . . , s K and let h be the pairwise distances between the s j 's. We shall denote by ν F λ (h) the empirical version of ν F λ (h), that is for s -

s p = h, ν F λ (h) = 1 2N N i=1 |F λ {Z i (s )} -F λ {Z i (s p )}|.
For a fixed, let θ a N LS (h) = ( θ a X (h), θ a Y (h)) t be estimated as above with some chosen distinct values λ ∈ Λ. Let λ / ∈ Λ and denote by H ⊂ [0, ∞) a finite set of spatial lags h. Moreover, let ν F λ be the estimation of ν F λ where θ X and θ Y are replaced in (3.1) by θ a X and θ a Y . We define:

(3.5) DC(a) = h∈H ω(h) ν F λ (h) ν F λ (h) -1 2
where the ω(h)'s are nonnegative weights which can be used for example to reduce the number of pairs included in the estimation. A simple choice for these weights is ω(h) = 1 {h≤r} where the r value can be chosen as the q-quantile of the distributions of the distances h between pairs of sites, q ∈ (0, 1). Finally, we select a that gives the lower value of (DC).

Of course, when dealing with real data, the marginal laws are usually not unit Frchet and thus have to be changed to unit Frchet. In that case, the empirical distribution function F is used instead of F in the definitions of Q i (h, λ) and ν λ F :

Q i (h, λ) = 1 2 | F λ (Z i (s)) -F λ (Z i (s + h))|, and
ν F λ (h) = 1 2N N i=1 | F λ {Z i (s )} -F λ {Z i (s p )}|.
We shall denote θ a X , θ a Y , a the estimations of θ X , θ Y and a obtained when Q i (h, λ) and ν F λ are used.

In order to get the consistency of our estimations, we need the two following assumptions:

• I 1 : for any a ∈ [0, 1],{λ 1 , λ 2 } ∈ Λ with λ 1 = λ 2 , the mapping [1, 2] 2 -→ R 2 (x, y) → (Φ(a, λ 1 , x, y), Φ(a, λ 2 , x, y))
is injective.

• I 2 : let θ X , (resp. θ Y ), θ X , (resp. θ Y ) be extremal coefficients of max-stable (resp. inverse max-stable) processes. Let λ be fixed, if for all h, Φ(a, λ, θ

X (h), θ Y (h)) = Φ(a , λ, θ X (h), θ Y (h)) then a = a .
Remark 3.2. The hypothesis I 1 and I 2 are identifiability hypothesis. Numerical tests on several models seem to indicate that they are satisfied for various max-mixtures models but we did not succeed to prove it. Theorem 3.3. Assume that (Z i (s j )) i=1,...,N are i.i.d copies of Z(s j ), j = 1, . . . , K where Z is a max-mixture spatial process with mixing coefficient a 0 ∈ [0, 1]. Assume that the injectivity conditions I 1 and I 2 are verified. Then the estimations of a by a and a are consistent in the sense that a -→ a 0 in probability as N → ∞ and a -→ a 0 in probability as N → ∞.

Proof. We shall give the proof for a, the proof for a is simpler and can be done along the same lines. We first begin by proving the consistency of θ a 0 X and θ a 0 Y . Consider {λ 1 , λ 1 } ∈ Λ, with

λ 1 = λ 1 . Write ε 1 h,i = Q i (h, λ 1 ) -Φ(a 0 , λ 1 , θ X (h), θ Y (h)) and ε 2 h,i = Q i (h, λ 1 ) -Φ(a 0 , λ 1 , θ X (h), θ Y (h)).
Using the convergence results from [START_REF] Naveau | Modelling pairwise dependence of maxima in space[END_REF]) (Proposition 3), we have that

1 N N i=1 ε k h,i → 0 in probability, k = 1, 2 and ∃ σ k ≥ 0 such that 1 N N i=1 (ε k h,i ) 2 → σ 2 k , in probability, k = 1, 2.
From this remark, following the lines of proof of Theorem II.5.1 in [START_REF] Antoniadis | Régression non linéaire et applications[END_REF], see also the proof of Theorem 4.1 in [START_REF] Ahmed | A semiparametric estimation for max-mixture spatial processes[END_REF], we conclude, using the injectivity hypothesis I 1 that for any h, θ a 0 X (h) -→ θ X (h) in probability and

θ a 0 Y (h) -→ θ Y (h) in probability. Now, consider (a * , θ * X (h), θ * Y (h)), a limit point of ( a, θ a X (h), θ a Y (h))
. Since a reaches the minimum of DC, we have DC( a) ≤ DC(a 0 ). The convergence of θ a 0 X , θ a 0 Y and ν F λ from [START_REF] Naveau | Modelling pairwise dependence of maxima in space[END_REF]

(Proposition 3) implies that lim N →∞ DC(a 0 ) = 0 so that DC(a * ) = 0 which leads to Φ(a * , λ , θ * X (h), θ * Y (h)) = Φ(a 0 , λ 2 , θ X (h), θ Y (h))
with w(h) = 0 and thus a * = a 0 by the injectivity condition I 2 .

Remark 3.4. As seen by our selection scheme above, we have to determine the choices of λ ∈ Λ and λ / ∈ Λ. After a lot of trials and as a compromise between accuracy and computation time we found that good results can be acheived with Λ = {1, 3} and λ = 1.5 as shown by the simulation study, Section 4.

Simulation Study

This section is devoted to a simulation study, in order to benchmark our estimation procedure. We shall consider two different max-mixture models.

M 1 is a max-mixture model in which X is a TEG process with A X a disk of fixed radius r X .

The AI process Y is an inverse TEG process with A Y a disk of fixed radius r Y . For simplicity, we choose stationary and isotropic exponential correlation functions, with range parameters

φ X , φ Y > 0 respectively. The model parameters vector is ψ = (r X , φ X , a, r Y , φ Y ) t .
M 2 is a max-mixture model where X is an isotropic Brown-Resnick process with variogram 2γ(h) = (h/φ X ) τ and Y is an isotropic inverted extremal-t process with v degrees of freedom and exponential correlation function

ρ Y (h) = exp(-h/φ Y ), φ Y > 0. The model parameters vector is ψ = (φ X , τ, a, v, φ Y ) t .
We summarize our simulation study procedure in the following steps.

Step 1. For each experiment, we consider a moderately sized data from the two max-mixture models described above with a true mixing coefficient a 0 , K = 50 sites randomly and uniformly distributed in the square [0, L] 2 , L ∈ N and N = 2000 independent replications at each site.

Max-stable processes were simulated using SpatialExtremes package in R [START_REF] Ribatet | Spatialextremes: modelling spatial extremes. R package version[END_REF] except for the TEG processes which have been simulated as in [START_REF] Anthony | [END_REF]. Each experiment was repeated M = 100 times.

Step 2. For each data set in Step 1, we estimated the extremal dependence functions θ X (h) and θ Y (h) using the nonlinear least squares estimation criterion as described in (3.4), using Q i .

This step was performed with a set of different mixing coefficients including the true one a 0 .

Step 3. We calculate our decision criterion (DC(a)) with the estimated θ a X (h) and θ a Y (h) functions from Step 2. It is expected that this criterion will lead to the minimum values when the true parameter a is used to estimate θ X (h) and θ Y (h) in Step 2. Since choosing a wrong mixing coefficient a should lead to a bad estimation of θ X (h) and θ Y (h). Equal weights ω(.) are used.

To assess the performance of the nonlinear least squares estimator θ a N LS (h), we simulate data from max-mixture models as mentioned in Step 1. The performances of the estimators are given by the relative mean square error MSE rel , for the K(K-1) 2 pairwise distances h j which are the distances between sites pairs (s , s k ): see [START_REF] Bel | Assessing extremal dependence of environmental spatial fields[END_REF] p. 171 for a similar definition in the multivariate context

(4.1) MSE rel (h j ) = M -1 M i=1 ( θ a (h j ) -θ(h j )) 2 θ(h j ) , j = 1, ..., K(K -1) 2 .
Figures 7 -10 display the estimation performances with respect to distance h using the mixing coefficients a = 0.75, 0.5, 0.25. For other max-mixture models we have similar results which are presented in the Supplementary Material (see Appendix B). Generally, our estimation procedure appears to work well for all considered max-mixture models. Moreover, the estimation of θ X (h) becomes more precise as the mixing coefficient value increases, and vice versa for θ Y (h). . Estimation performance of θNLS (h) estimates. Data simulated from max-mixture model M 1 with parameter ψ = (0.2, 0.1, 0.25, 0.9, 0.7) t in the square [0, 2] 2 as described in Step 1. The top row compares the true extermal coefficient functions (black points) and the NLS mean estimates (blue points) against the distance h. The middle row displays the MSE rel of the NLS estimates against the distance h. While the bottom row displays the histograms of MSE rel of the NLS estimates.

0 1 2 3 4 5 6 1.0 1.4 1.8 True θ X (h), mean of θ X (h) h 0 1 2 3 4 5 6 1.0 1.4 1.8 True θ Y (h), mean of θ Y (h) h 0 1 2 3 4

Real data example

The data analysed in this section are daily rainfall amounts in (millimetres) over the years Following the approach of [START_REF] Bacro | A flexible dependence model for spatial extremes[END_REF] graphical assessments to explore possible anisotropy Figure 13 based on the empirical estimates of the functions χ(h, u) and χ(h, u) We apply our methodology for the selection of the mixing coefficient for all a ∈ (0, 1) by steps 0.01. The a → DC(a) function is plotted in Figure 14. The best-fitting max-mixture model as judged by our DC criterium has a mixing coefficient a = 0.34. In the literature, the widely used parametric inference procedure is based on composite likelihood methods. In particular, pairwise likelihood estimation has been found usefull to estimate parameters in a max-stable process. A description of this method can be found in [START_REF] Simone A Padoan | Likelihood-based inference for maxstable processes[END_REF][START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF][START_REF] Bacro | A flexible dependence model for spatial extremes[END_REF] for spatial context.

Unfortunately, parameter estimation using composite likelihood suffers from some defects.

First, it can be onerous, since the computation and subsquent optimization of the objective function is time-consuming. Second, the choice of good initial values for optimization of the composite likelihood is essential. Third, model dependency, a preliminary step to conduct a composite likelihood estimation is to specify the model that describes the dependence structure. So, this mission seems to be laborious, due to the large number of combinations that can be formed from the AD and AI processes stemming from max-stable processes, since the variety of dependence structures that can be assumed, i.e. changing the correlation coefficient function type in Schlather model [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] leads to different dependence structures. So, an unacurate choice may lead to severe under/over estimations of probabilities associated to simultaneous extreme events.

In the sequel, we shall estimate the conditional probability of having daily rainfall that exceeds some threshold z at an unused site denoted by s * 0 given that this event has occurred at the nearest observed site which is denoted by s 0 , i.e., P[Z(s * 0 ) > z|Z(s 0 ) > z]. We compare this estimation with that obtained by the best-fitting parametric model based on composite likelihood estimation.

For this purpose, we fitted the daily rainfall data based on censored pairwise likelihood approach used by [START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF][START_REF] Bacro | A flexible dependence model for spatial extremes[END_REF] where the threshold is taken corresponding to the 0.9 empirical quantile at each site. We fitted the generalized extreme value distribution GEV(µ, σ, ξ) separately to each site and then data are transformed to unit Fréchet margins through the transformation z → -1 log( G(z))

, where G(.) is the estimated GEV cumulative distribution function. The models are M a : a MM model where X is a TEG process with an exponential correlation function ρ(h) = exp(-h /φ X ), φ X > 0. A X is a disk of fixed and unknown radius r X , and Y is an inverted TEG process with exponential correlation function ρ(h) = exp(-h /φ Y ), φ Y > 0, and A Y is a disk with fixed and unknown radius r Y .

M b : a MM model where X is a TEG process as in M a . Y is an isotropic inverted Smith process where Σ is a diagonal matrix (σ 12 = 0) with σ 2 11 = σ 2 22 = φ 2 Y , i.e., γ(h) = ( h /φ Y ).

M c : a MS TEG process described as X in M a .

M d : a MS isotropic Smith process where Σ is a diagonal matrix (σ 12 = 0) with σ 2 11 = σ 2 22 = φ 2 X , i.e., γ(h) = ( h /φ X ).

Conclusion

In this paper, we have proposed a statistically efficient nonparametric model-free selection criterion. We can exploit our decision about the mixing coefficient to predict conditional probabilities of daily rainfall at unobserved sites depending on the dependence structure in the analyzed data.

The main advantage of our approach is that it is model free, unlike the parametric approach which assumes a specified model, so the risk of unaccurate choice of stochastic processes for describing the joint tail distribution may lead to severe under/over estimation of probabilities associated to simultaneous extreme events.

We have shown in our real data example that the max-mixture approach appears of interest for modeling environmental data. In particular it has the eligibility to overcome the limitations of the max-stable models in which only asymptotic dependence or exact independence can be modeled. 2.12) on the logarithm scale according different values of mixing coefficient a ∈ {1, 0.75, 0.5, 0.25, 0}. Top row: X is isotropic extremal-t with v 1 = 1 degrees of freedom and ρ(h) = exp(-h), Y is isotropic inverted extremal-t process with v 2 = 2 degrees of freedom and ρ(h) = exp(-h/1.5). bottom row: X is isotropic Brown-Resnick with variogram 2γ(h) = h 2 , X is isotropic inverted Brown-Resnick with variogram 2γ(h) = (h/1.5) 2 .

  where ε k (s) are independent replicates of a stationary Gaussian process ε = {ε(s), s ∈ S} with zero mean, unit variance and correlation function ρ(.). 1 A is the indicator function of a compact random set A ⊂ S, A k are independent replicates of A and R k are points of a Poisson process with a unit rate on S. The constant c is chosen to satisfy the constraint E{Q k (s)} = 1.

Figure 1 .

 1 Figure1. Simulations of two inverted max-stable processes (2.11) on the logarithm scale. Left panel: simulated image of AI process constructed by inverting an isotropic inverted Brown-Resnick process with variogram 2γ(h) = (h/1.5) 2 . Middle panel: simulated image of AI process constructed by inverting an isotropic extremal-t process with v = 1 degrees of freedom and exponential correlation function ρ(h) = exp(-h/1.5). On the Right panel: associated functions η(h).

Figure 2 Figure 2 .

 22 Figure 2 displays five simulated images of the max-mixture model over the [0, 10] 2 square according to different values of the mixing coefficient a. In order to show the role of the mixing

Figure 3

 3 Figure 3 illustrates how the spatial extremal dependencies vary regarding to a. Based on

Figure 3 .

 3 Figure3. Associated θ(.) functions for a max-mixture process Z(.) based on two max-stable processes X(.) and X (.) that belong to the same family with the same parameters used in Figure1for a ∈ {0.25, 0.5, 0.75}.

Figure 4 Figure 4 .

 44 Figure4displays the behavior of the F λ -madogram for three different max-mixture models with respect to the distance h for different λ ∈ {0.5, 1, 1.5, 3}. λ = 1.5 corresponds to the

Figure 5 .Figure 6 .

 56 Figure5. Estimation performance of θNLS (h) estimates. Data simulated from max-mixture model M 1 with parameter ψ = (0.2, 0.1, 0.75, 0.9, 0.7) t in the square [0, 2] 2 as described in Step 1. The top row compares the true extermal coefficient functions (black points) and the NLS mean estimates (blue points) against the distance h. The middle row displays the MSE rel of the NLS estimates against the distance h. While the bottom row displays the histograms of MSE rel of the NLS estimates.

Figure 9 .

 9 Figure 9. Estimation performance of θNLS (h) estimates. Data simulated from max-mixture model2 with parameter ψ = (0.1, 2, 0.5, 2, 1.5) t in the square [0, 5] 2 as described in Step 1. The top row compares the true extermal coefficient functions (black points) and the NLS mean estimates (blue points) against the distance h. The middle row displays the MSE rel of the NLS estimates against the distance h. While the bottom row displays the histograms of MSE rel of the NLS estimates.

Figure 11 .

 11 Figure11. Top row: boxplots for decision criterion of 100 data replications of 2000 independent copies from M 1 with parameters φ X = 0.1, r X = 0.2, φ Y = 0.7 and r Y = 0.9 over the square [0, 2] 2 . Bottom row: boxplots for decision criterion of 100 data replications of 2000 independent copies from M 2 with parameters φ X = 0.1, τ = 2, v = 2 and φ Y = 1.5 over the square [0, 5] 2 .

  1972-2016 occurring during April-September at 38 sites in the East of Australia whose locations are shown in Figure 12. The altitude of the sites varying from 4 to 552 meters above mean sea level. The sites are separated by distances of approximately (34-1383)km. This data set is freely available on http://www.bom.gov.au/climate/data/index.shtml?bookmark=136.

Figure 12 .

 12 Figure 12. Geographical locations of 41 meteorological stations located in the East of Australia. Black crosses represent the 38 stations used for selecting the mixing coefficient a. The three colored numbers {1,2,3} correspond to the unused stations on which we shall predict conditional probabilities in order to validate our procedure.

Figure 14 .

 14 Figure14. Decision criterion values for real data example on the interval a ∈ (0, 1) by steps 0.01. The blue star corresponds to the minimum DC value which occurs at a = 0.34.

Figure 16 .

 16 Figure16. Simulations of the max-mixture model (2.12) on the logarithm scale according different values of mixing coefficient a ∈ {1, 0.75, 0.5, 0.25, 0}. Top row: X is isotropic extremal-t with v 1 = 1 degrees of freedom and ρ(h) = exp(-h), Y is isotropic inverted extremal-t process with v 2 = 2 degrees of freedom and ρ(h) = exp(-h/1.5). bottom row: X is isotropic Brown-Resnick with variogram 2γ(h) = h 2 , X is isotropic inverted Brown-Resnick with variogram 2γ(h) = (h/1.5) 2 .

Figure 17 .

 17 Figure17. Theoretical F λ -madogram functions 3.1 with λ ∈ {0.5, 1, 1.5, 3}. Left panel: Max-mixture model in which in which X is isotropic Brown-Resnick with variogram 2γ(h) = (h/0.5) 2 , and Y is isotropic inverted extremal-t process with v = 1 degrees of freedom and ρ(h) = exp(-h/1.5). Middle panel: Max-mixture model in which X is isotropic extremal-t with v 1 = 1 degrees of freedom and ρ(h) = exp(-h), Y is isotropic inverted extremal-t process with v 2 = 2 degrees of freedom and ρ(h) = exp(-h/1.5). Right panel: Max-mixture model in which X is isotropic Brown-Resnick with variogram 2γ(h) = (h/0.2) 2 , Y is isotropic inverted Brown-Resnick with variogram 2γ(h) = h 2 .

  Figure8. Estimation performance of θNLS (h) estimates. Data simulated from max-mixture model M 2 with parameter ψ = (0.1, 2, 0.75, 2, 1.5) t in the square [0, 5] 2 as described in Step 1. The top row compares the true extermal coefficient functions (black points) and the NLS mean estimates (blue points) against the distance h. The middle row displays the MSE rel of the NLS estimates against the distance h. While the bottom row displays the histograms of MSE rel of the NLS estimates.

	MSErel(θX(h))	0.000 0.008						MSErel(θY(h))	0.01 0.04
							5	6	0	1	2	3	4	5	6
					h							h
					MSE rel (θ X (h))					MSE rel (θ Y (h))
	Frequency	0 400 800						Frequency	0 200 400
		0.000	0.002	0.004	0.006	0.008	0.010	0.012	0.00	0.01	0.02	0.03	0.04	0.05	0.06
					MSErel							MSErel
				True θ X (h), mean of θ X (h)				True θ Y (h), mean of θ Y (h)
		1.0 1.4 1.8							1.0 1.4 1.8
		0	1	2	3	4	5	6	0	1	2	3	4	5	6
					h							h
	MSErel(θX(h))	0.00 0.02 0.04						MSErel(θY(h))	0.005 0.020
		0	1	2	3	4	5	6	0	1	2	3	4	5	6
					h							h
					MSE rel (θ X (h))					MSE rel (θ Y (h))
	Frequency	0 400 1000						Frequency	0 150 300
		0.00		0.01	0.02		0.03	0.04	0.000	0.005	0.010	0.015	0.020	0.025	0.030
					MSErel							MSErel

  Figure10. Estimation performance of θNLS (h) estimates. Data simulated from max-mixture model M 2 with parameter ψ = (0.1, 2, 0.25, 2, 1.5) t in the square [0, 5] 2 as described in Step 1. The top row compares the true extermal coefficient functions (black points) and the NLS mean estimates (blue points) against the distance h. The middle row displays the MSE rel of the NLS against the distance h. While the bottom row displays the histograms of MSE rel of the NLS estimates.

			True θ X (h), mean of θ X (h)					True θ Y (h), mean of θ Y (h)
		1.0 1.4 1.8							1.0 1.4 1.8
		0	1	2	3	4	5	6	0	1	2	3	4	5	6
					h							h
	MSErel(θX(h))	0.00 0.04						MSErel(θY(h))	0.000 0.006
		0	1	2	3	4	5	6	0	1	2	3	4	5	6
					h							h
					MSE rel (θ X (h))						MSE rel (θ Y (h))
	Frequency	0 100 250						Frequency	0 100 250
		0.00	0.02		0.04	0.06		0.08	0.000	0.002	0.004	0.006	0.008	0.010
					MSErel							MSErel

M e : the inverted Smith process described as Y in M b .

The composite likelihood information criterion (CLIC) [START_REF] Varin | A note on composite likelihood inference and model selection[END_REF], defind as CLIC= -2[p ( ϑ)-tr{J ( ϑ)H -1 ( ϑ)}] is used to judge the best-fitting model. Here, the maximum pairwise likelihood estimator is denoted by ϑ. Lower values of CLIC indicate a better fit. Our results are summarised in Table 1. The best-fitting model for the data, as judged by CLIC, is the hybrid dependence model M b .

where h 0 is the separation distance between s * 0 and s 0 . Equation (5.1) may be used to estimate

using both parametric and nonparametric approaches. We consider s * 0 as the three unused sites that have been marked by colored numbers {1,2,3} on the map Figure 12.

The threshold z is taken corresponding to the q-emperical quantile at the site s 0 , q ∈ (0, 1).

For estimating P[Z(s * 0 ) > z|Z(s 0 ) > z] nonparametrically, we fitted the data again by NLS procedure with the best mixing coefficient a = 0.34 and we obtained the estimators θ X (h 0 ) and θ Y (h 0 ) by averaging the values of θ X (h) and θ Y (h) , where h ∈ [h 0 -10, h 0 +10]km by taking the advantage of stationarity and isotropy of our data and h 0 ∈ {54.226, 92.534, 133.673}km. While, for estimating P[Z(s * 0 ) > z|Z(s 0 ) > z] parametrically, θ X (h) and θ Y (h) are obtained by substituting the separating distance h 0 and the estimated parameters from hybrid model M b Table 1.

In order to compare the results obtained by the two approaches, we used the data at the unused three sites to compute the empirical version of conditional probabilities P[Z(s * 0 ) > z|Z(s 0 ) > z]. Below are the resulting P-P plots. Figure 15. Diagnostic P-P plots for threshold excess conditional probabilities for the three unused sites {1,2,3} on the map Figure 12 obtained by both approaches. The best parametric model M b as judged by the CLIC and our nonparametric approach. Green: site 1; red: site 2; blue: site 3.

Generally, Figure 15 shows that our nonparametric approach outperforms the parametric one for predicting P[Z(s * 0 ) > z|Z(s 0 ) > z]. One of the justifications for this situation is that with the parametric model we have to specify a model that describes the dependence structure, and listing all choices seems to be a tedious task due to the large number of possibilities that can emerge from the AD and AI processes stemming from max-stable processes, i.e., different choices of correlation coefficient functions for the same model leads to different models. So, in this case the inaccurate choice/guess of models to be fitted may lead to severe under/over estimation of probabilities associated to simultaneous extreme events. Supplementary Material related to this article can be found online at http://math.univ-lyon1.fr/homes-www/abuawwad/Supplementary/.

Appendix C.

Proof of Lemma 5.1. Denoting the joint survivor function of the process Z (2.12) P[Z(s 1 ) >

x, Z(s 2 ) > y] by ḠZ (x, y), we may write

where G X (x), G Y (y) are the marginal probability distribution functions of processes X(.) and Y (.) respectively in model (2.12) and G Z (x, y) is the bivariate probability distribution function of the stochastic process Z(.) in the same model. using (2.14) and (2.16), setting x = y = z

where ω 1-a z = -1/ log 1 -e -1-a z , taking the advantage of the homogeneity of order -1 of the exponent measures V X (.) and V Y (.); we have V X (z, z) = θ X /z, e [-V Y [ω( 1-a z ),ω( 1-a z )] =

1 -e -1-a z θ Y and this gives (5.1).


Benoît Vandame
email: [benoit.vandame@technicolor.com

Valter Drazic
email: valter.drazic@technicolor.com

Matthieu Hog
email: matthieu.hog@technicolor.com

Neus Sabater
email: neus.sabater]@technicolor.com

Plenoptic sensor: application to extend field-of-view

In this paper we study the light field sampling produced by ideal plenoptic sensors, an emerging technology providing new optical capabilities. In particular, we leverage its potential with a new optical design that couples a pyramid lens with an ideal plenoptic sensor. The main advantage is that it extends the field-of-view (FOV) of a main-lens without changing its focal length. To evince the utility of the proposed design we have performed different experiments. First, we demonstrate on simulated synthetic images, considering ideal and real lenses, that our optical design effectively doubles the FOV. Then, we show its feasibility with two different prototypes using plenoptic cameras on the market with very different plenoptic samplings, namely a Raytrix R5 and a Canon 5D MarkIV. Arguably, future cameras with ideal plenoptic sensors will be able to be coupled with pyramid lenses to extend its inherent FOV in a single snapshot.

Introduction

Plenoptic cameras are able to capture the Light-Field (LF), thanks to a micro-lens array (MLA) placed between the main-lens and the sensor. Depending on the MLA position, plenoptic cameras are divided in type-1 [START_REF] Ng | Digital Light Field Photography[END_REF] (e.g. Lytro [1]) and type-2 or focused [START_REF] Lumsdaine | The focused plenoptic camera[END_REF] (e.g. Raytrix [3]).

The first plenoptic cameras on the market have the particularity that the MLA is not well aligned with the sensor. Indeed, plenoptic camera manufacturers assemble individual optic components producing an unavoidable rotation offset between the MLA and the pixel matrix. For this reason, many research works in the field have focused on calibration and decoding methods [START_REF] Dansereau | Decoding, calibration and rectification for lenselet-based plenoptic cameras[END_REF][START_REF] Hog | An image rendering pipeline for focused plenoptic cameras[END_REF], as well as image processing algorithms taking into account such misalignments [START_REF] Cho | Modeling the calibration pipeline of the lytro camera for high quality lightfield image reconstruction[END_REF][START_REF] Sabater | Accurate disparity estimation for plenoptic images[END_REF][START_REF] Fiss | Refocusing plenoptic images using depth-adaptive splatting[END_REF]. However, camera manufacturers rely on wafer-level fabrication to assemble micro-optical components like MLA's onto pixel matrix with great accuracy. Thus, plenoptic cameras with the MLA ideally aligned with respect to the sensor shall be available in the near future Such ideal plenoptic sensors (simply called plenoptic sensors in the sequel) provide a new paradigm for LF processing since no camera calibration is required and sub-aperture images (SAI's) or Epipolar Plan Images (EPI's) are simply extracted without interpolation.

In fact, very simple plenoptic sensors already exist in the consumer market. This is the case of dual-pixels in highend smartphones such as the Samsung S7 [START_REF] Choi | An all pixel pdaf cmos image sensor with 0.64µmx1.28µm photodiode separated by self-aligned in-pixel deep trench isolation for high af performance[END_REF], and DSLR (Digital Single Lens Reflex) cameras such as the Canon 5D MarkIV. These devices provide a limited angular sampling of the LF comparable to the right and left views of a stereo camera. Maintaining a zero offset between these two views controls the autofocus of the camera.

Wide-FOV imaging is achieved by stitching multiple images that are recorded from the same center of projection [START_REF] Szeliski | Creating full view panoramic image mosaics and environment maps[END_REF]. Stitched images produce a better spatial resolution, but the parallax between the views produce artifacts on the resulting image. This issue is addressed by the LF panorama stitching [START_REF] Zhou Xue | Light field panorama by a plenoptic camera[END_REF]. Unfortunately, both strategies require a sequential capture, thus dedicated to static scenes. Changing the main-lens for a wider FOV lens is another solution, but these lenses are often bulky and require larger stack of lenses to produce sharp images. Alternatively, monocentric lenses [START_REF] Stamenov | Optimization of two-glass monocentric lenses for compact panoramic imagers: general aberration analysis and specific designs[END_REF][START_REF] Stamenov | Panoramic monocentric imaging using fiber-coupled focal planes[END_REF] have become increasingly popular for gigapixel imaging [START_REF] Cossairt | Gigapixel computational imaging[END_REF]. With these lenses, the light is collected on a spherical surface either with a curved sensor or a fiber coupling interface to a flat sensor. This approach has been pushed further using a plenoptic camera [START_REF] Dansereau | A wide-field-of-view monocentric light field camera[END_REF]. Finally, combining a prism array with common cameras has been proposed to double [START_REF] Chen | Optics system design applying a micro-prism array of a single lens stereo image pair[END_REF] or quadruple [START_REF] Sun | Singlelens camera based on a pyramid prism array to capture four images[END_REF] the FOV of the main-lens. In [START_REF] Somanath | Single Camera Stereo System Using Prism and Mirrors[END_REF], prisms and mirrors are combined for stereo capture out of a single lens camera. Our contributions In this paper we analyze the advantages and constraints of plenoptic sensors. In particular, we focus on the so-called quad-pixel sensor where a micro-lens covers 2 × 2 pixels. In our study, we describe the relationship between the SAI's and the corresponding portions of the main-lens pupil through which light rays have travelled.

Besides, we propose an optical design for a single lens camera that doubles the FOV of the lens combining a plenoptic sensor and a pyramidal lens made of four prisms. The main idea is that the prisms deviate the photons entering the main-lens creating four distinct views that can be demultiplexed into the SAI's thanks to the plenoptic sensor. Then, with a single snapshot, the stitching of the SAI's increases the FOV up to a factor of two in each direction without changing the main-lens focal length, which is an unprecedented capability. Furthermore, the use of plenoptic sensors simplifies the parameterization and the processing of the captured LF compared to existing designs in the literature. Our experiments include synthetic image simulations and real images captured with two different prototypes we have constructed.

From plenoptic cameras to plenoptic sensor

In this paper we focus on type-1 plenoptic sensors for its capacity to sample the main-lens exit pupil. Indeed, type-1 plenoptic cameras [START_REF] Ng | Digital Light Field Photography[END_REF] are characterized by the fact that the distance d between the MLA and the sensor is equal to the micro-lenses focal length f (as illustrated in Fig. 1). In that case, considering a thin lens model, the micro-lenses focus at infinity. Now, the main-lens is considered at the microlenses optical infinity since (D d). So, the micro-lenses are imaging the main-lens exit pupil which is equal to the lens aperture in a thin lens model. The image on the pixels underneath the micro-lens is called micro-image.

Captured micro-images with N × N pixels, N ∈ N, correspond to sampling the main-lens exit pupil with a N × N grid. For instance, if N = 2, each pixel of the 2 × 2 micro-image integrates the light rays passing through one of the four quarter-discs of the aperture stop.

It turns out that the sharp sampling of the exit pupil with the same pixel grid for all micro-images requires the MLA to be perfectly aligned with the sensor, this is the same squared lattice and same orientation than the pixel matrix. By contrast, in type-2 or focused plenoptic cameras (d = f ) micro-lenses focus at a plane that does not match with the main-lens exit pupil [START_REF] Georgiev | Focused plenoptic camera and rendering[END_REF].

Regarding the optical design of type-1 plenoptic cameras, it is imposed that the F-numbers of the main-lens and micro-lenses are equal. In that case, the system is said to be aperture matched and it guarantees the micro-images to cover as many pixels as possible without overlapping. Indeed, using the Thales theorem, the distance P between two consecutive micro-images (in physical unit) is equal to P = φ e, where e = 1 + d/D. Tuning the main-lens Fnumber, F/Φ, the diameter P of the micro-images is set to be equal to P ≈ φ.

Pupil sampling with the micro-images

Ideally the distance P between two micro-images should be equal to an integer number of pixels. This is, P = N δ, with δ being the pixel size. This requires the micro-lenses to have a pitch φ = N δ/e which is slightly smaller than N pixels and is function of the main-lens distance D. Fig. 2 illustrates the micro-images position with respect to the pixel array at the center of the sensor and at the border.

It is worth mentioning that sensor manufacturers use micro-lenses to guide photons to the middle of each pixel where the photo-diode is located. Recently, manufacturers have also designed a matrix of pixels mounted with a MLA having a pitch φ slightly smaller than the pixel size δ [START_REF] Sekine | Solid state image sensor, method of manufacturing solid state image sensor, and image capturing system[END_REF]. It corresponds to the case N = 1 with all the micro-images centered on the middle of the corresponding pixels. In other words, the chief rays, i.e. the rays passing through the exit pupil and the micro-lens centers, hit the photo-diodes despite the increasing Chief Ray Angle (CRA), which is the angle between a chief ray and the optical axis of a microlens. The so-called CRA correction decreases the pixel vignetting [START_REF] Peter Bert Catrysse | Qe reduction due to pixel vignetting in cmos image sensors[END_REF]. This fact emphasizes how ready are manufacturers to produce quad-pixels plenoptic sensor with CRA correction.

Study of SAI's on the Plenoptic Sensor

Collecting the SAI's is simple thanks to the integer size N × N of the micro-images. Let L(x, y) be the image cap-

S i,j (k, l) = L x N + i mod N, y N + j mod N (1)
where

(k, l) = ( x/N , y/N ) ∈ [0, N x /N ] × [0, N y /N ].
If the micro-lens pitch φ is not exactly equal to N δ/e, then P is not an integer number of pixels, and computing the SAI's with Eq. 1 is incorrect. Indeed, by definition SAI's require to collect pixels at a fixed distance from the micro-image centers. So, if P = N δ, SAI computation requires interpolation to extract pixels at non-integer coordinates from the sensor image [START_REF] Dansereau | Decoding, calibration and rectification for lenselet-based plenoptic cameras[END_REF]. Such interpolation averages micro-image pixels, mixing the angular information encoded in the micro-images. To prevent it, Eq. 1 is nevertheless used to extract approximate SAI's that we note S i,j , even if P is not a multiple of δ. Such approximation amounts to interpolate with nearest neighbors instead of a more sophisticated method.

It is interesting to point out that the micro-images sample the main-lens exit pupil with a constant partition when P = N δ regardless of the micro-image position on the sensor. On the contrary, when P = N δ, the approximate SAI's do not sample the main-lens exit pupil homogeneously. Fig. 3 illustrates the portion of the pupil visible by the SAI S (0,0) of a quad-pixel (N = 2) type-1 plenoptic sensor.

Considering real-lenses

Until now, we have modelled the main-lens with a thin lens model. However, considering the real-lens is mandatory to understand how the main-lens exit pupil is sampled. The real-lens pupil, also named the aperture stop, is the physical stop delimiting the beam of photons entering the camera. The aperture stop is located where a diaphragm can reduce homogeneously the amount of light on the sensor. Most lenses designed for DSLR cameras have an aperture stop roughly located within the main-lens (as illustrated by a typical double-Gauss lens in Fig. 4-left). On the contrary, for smartphones, the aperture stop is mostly located at the entrance of the main-lens [START_REF] Mercado | Small form factor telephoto camera[END_REF] (see Fig. 4-right), even though smartphones do not have diaphragms. In this paper, we consider main-lenses where the aperture stop is positioned at the first diopter on the opposite side of the sensor.

With real-lenses, the distance D is the distance between the so called exit pupil and the MLA. The position of the exit pupil is located at the imaging plane of the aperture stop by the rear lenses and depends on the individual lenses located between the aperture stop and the sensor.

Doubling the FOV

Extending the FOV with a pyramid lens With a common camera, θ the half FOV of an image is given by the focal length F of the main-lens and the physical size T of the sensor: θ = arctan T 2F . To extend the FOV, a pyramid lens (i.e. four prisms assembled together) is placed at the main-lens aperture stop. The FOV per portion of the pupil is rotated in different orientations.

Each prism deviates the photons entering the main-lens by an angle α with respect to the x axis. α is chosen to be equal to the half FOV angle θ of the main-lens. With the two prisms of Fig. 5, the sensor records the superposition of two images, each one imaging a different part of the mainlens exit pupil. These two superposed images, combined together, double the FOV in the vertical direction of the camera. Similarly, using four prisms in front of the main-lens, each one covering a quarter of the main-lens aperture stop, doubles the FOV of the camera in vertical and horizontal orientations.

The angle α associated with the prism is function of its angle A and the refraction index n of its material. A good approximation gives α ≈ (n -1)A. Since α is set to be equal than θ, the prism angle A is easily computed knowing the refraction index of the prism material.

Extending the FOV with a plenoptic sensor Discriminating the two or four images that have been summed at the sensor plane is not easy. It requires to know which prism the photon has crossed. This is done with a plenoptic sensor with an ideal pupil sampling. Considering a quad-pixel sensor, the four SAI's are stitched to produce a single image that has a double FOV, as many pixels as the sensor and a single exposure time. Note that the equivalent f-number of the stitched image is divided by two, since SAI's collect only a quarter of the incoming photons.

Image simulation

Simulated images are generated with PBRT [START_REF] Pharr | Physically Based Rendering, Second Edition: From Theory To Implementation[END_REF] which has been extended to support real lenses. Our extension permits to define: a thin lens model, real-lenses with their interfaces (each one is defined with a mathematical model of its curvature and a glass material similar to [START_REF] Kolb | A realistic camera model for computer graphics[END_REF]), an array of prisms and a MLA. In particular, a ray is cast within a cone which is defined by an apex located on the sensor and a disk-shaped base which is defined by the exit pupil diameter of the main-lens and is located at the exit pupil position from the sensor. PBRT defines 5D random coordinates: 2D for the sensor coordinates (cone apex), 2D for the exit pupil coordinates (within the cone apex) and 1D for the wavelengths. Using the Snell-Descartes law, a ray is refracted at each interface along its path.

Ideal Lens simulation

Primarily, synthetic images are generated with the ideal thin lens model (for both main-and micro-lenses). The plenoptic camera characteristics are defined in Table 1. The sensor size is T = 3.6 × 3.6mm 2 as typically found in smartphones. By design (P = φe), the size of the microimages is strictly equal to 2 × 2 pixels. A test chart of colored letters and numbers is located at z = 2m from the camera. The main-lens produces a sharp image on the micro-lenses, D is computed using the thin lens equation. Four prisms forming a pyramid are positioned at the mainlens aperture stop. Each prism is defined by a material with a constant refraction index of n = 1.74 and an angle of

F = 6.16mm Φ = 3.09mm D = 6.18mm f = 4.8µm φ = 2.3981µm d = f δ = 1.2µm N x N y = 3000 2 T = 3.6mm
Table 1: Ideal plenoptic camera characteristics. A = 16.7 • . Fig. 6 illustrates the SAI's extracted from the synthetic simulation. Without prisms the four SAI's are identical whereas with the four prisms, each SAI observes a deviated FOV. The observed distortions are due to the prism. Cancelling the constant geometrical distortion enables to stitch the four SAI's. The FOV of the main-lens was 32 • , and with the prisms becomes 64 • (or 77 • considering the hypotenuse of the sensor with a 2/3 ratio), competing with a typical wide angle lens for smartphones.

It is worth mentioning that the image sharpness of one SAI is quite equivalent to the image sharpness of the mainlens alone. The simulated main-lens is optimized for a field of 32 • , beyond that field, the image quality degrades dramatically. This experiment also demonstrates that narrowangle main-lenses with simple optical design can produce extended FOV with their native resolution.

Real-lens simulation

Using Zemax optics studio [START_REF]Optic studio[END_REF], we have designed and optimized a real main-lens made of three biconvex lenses. It has comparable characteristics as the ideal lens simulation shown above. The six interfaces (two interfaces per lens) are presented sequentially starting from the closest interface to the sensor (see where the optic axis is presumed to lie in the z direction, and z(r) is the sag (the z-component of the displacement of the surface from the vertex, at distance r from the axis). More precisely, the sag model is an aspheric surface characterized by a radius of curvature R and a conic term K:

z(r) = r 2 R 1 + (1 -K) r 2 R 2 . ( 2
)
An interface is fully defined by the thickness between itself to the next one measured on the main optical axis, a curvature radius (in mm), the physical material following the interface, a diameter (in mm) and a conic term.

The refractive index of a material is function of the wavelength. The function is often modeled by the Sellmeier equation made of 6 parameters [2]. In our simulation software, we use this model to estimate the refractive index of a given material for the three color channels RGB. Fig. 7 shows the six interfaces and the sensor. The interfaces have been computed such that the main-lens aperture stop is located at the first interface where the four prisms should be positioned. The geometrical distortion of this main-lens is within ±0.5%. The MLA and sensor are simulated with the same thin lens model than the ideal simulation. The four prisms are modelled with a perfect material with a constant refraction index (thus producing no chromatic aberrations).

The MLA is designed such that SAI's collect photons strictly from a corresponding quarter-pupil of the main-lens. However, with real-lenses, the rays are distorted due to the lens aberrations. As a result, some photons passing by a given quarter-pupil are not recorded by the expected SAI. We call this phenomenon sub-aperture cross-talk. To illustrate it, we simulate a white board located at z = 2m and a mask at the main-lens aperture stop to mask out three of the four quarter-pupils. One expects only one SAI to record photons but we have observed that this is not the case (see Fig. 8). The cross-talk ratio r is the ratio between the maximum flux on another SAI versus the flux on the corresponding SAI. We have measured that r = 13%, which is too high to completely isolate the photons passing from different prisms. To solve this problem, a cross shaped mask is positioned at the aperture stop of the main-lens. The thickness of the mask is set to have a ratio r ≈ 3%. The mask width is 1/15 of the aperture stop diameter, which removes 16% of the photon flux entering the main-lens. Note that the cross-talk ratio depends on the main-lens design. According to our simulations, it seems correlated with the geometrical distortion of the main-lens. Fig. 9 compares two SAI's of a test chart with and without the cross-shaped mask.

Extended FOV with a type-2 camera

We ideally would experiment with a type-1 camera, but the Raytrix camera is the only available plenoptic camera which enables to change the main-lens. Indeed, a main-lens with an aperture stop located at its entrance is required for our experiment. In particular, the experiment we have set up allows to double horizontally the FOV of a main-lens mounted with the Raytrix R5 camera.

Type-2 pupil sampling

As mentioned before, micro-lenses in type-2 cameras do not focus at the main-lens exit pupil but at a plane S distant by d from the MLA. Thus, the main-lens exit pupil is not sampled with a disjoint partition. Instead, corresponding pixel visibility areas overlap (see Fig. 10). The overlap thickness between two juxtaposed portions is function of the ratio f /d and the pupil diameter Φ.

Converting a type-2 into a dual-pixel camera

Dual-pixel refers to a plenoptic sensor having 2 rectangular pixels bellow each micro-lens. It enables to capture 2 SAI's. Our experiment converts the Raytrix camera into a dual-pixel camera. Our prototype has only one prism to deviate half of the photons entering the main-lens. The mainlens needs to be sufficiently thin, so the prism can be set as close as possible to it. We assume the prism to be placed at the main-lens aperture stop. Supplementary material illustrates the experiment with the Raytrix camera.

The Raytrix R5 has many pixels per micro-lens (P = 20.20 δ). To convert it into a dual-pixel camera, stereo refocused images are computed splitting the left and right portions of the micro-images. The two refocused images are comparable to the SAI's extracted from a dual-pixel sensor. We have used the pipeline in [START_REF] Hog | An image rendering pipeline for focused plenoptic cameras[END_REF] to compute the two refocused images from the Raytrix R5. We summarize the main steps bellow:

Micro-lens image calibration The MLA is an hexagonal lattice fully characterized by a radius, a rotation angle (with respect to the pixel matrix) and the offset between the first micro-image center and the origin of the sensor. These values are all computed using the Fourier transform of a white image, this is a shot of a flat white screen homogeneously illuminated. The calibration associates a pixel coordinate (x, y) with a micro-lens coordinate (i, j).

Micro-lens vignetting correction

The white image is normalized by its average maximum and inverted to define a flux scaled correction per pixel that fixes the vignetting of the micro-images.

Stereo image refocusing Image refocusing with a type-2 camera is defined with Eq. 3:

X Y = s g x y - x i,j y i,j + x i,j y i,j (3)
Where (X, Y ) is the projection on the refocused image of the sensor pixel (x, y) which belongs to micro-image (i, j). Micro-image (i, j) is centered at pixel coordinate (x i,j , y i,j ). g controls the refocalization distance, and s the relative size of the refocused image versus the input LF image. The LF pixel (x, y) is projected at the non-integer coordinate (X, Y ). Lanczos interpolation is used to splat the pixel into the refocused image. A weight-map is maintained to count the projected pixels. The refocused image is normalized with the weight-map after all pixels are projected.

The left and right refocused images are computed with Eq. 3 and the following conditions respectively (x-x i,j ) < -m/2 and (x -x i,j ) > m/2. Where m is the thickness in pixels of the vertical masks located at the middle of the micro-images to cancel sub-aperture cross-talk.

The refocalization parameters are set to s = 0.2 (refocused images are 5 times smaller than the input image) and g = 4. The test chart is positioned such that the images appear as sharp as possible. Several prisms with different angles are tested such that the left and right images are almost juxtaposed with a tiny common vertical portion. Experimentally m is set to 2 to cancel the cross-talk.

Extended image stitching The prism produces small distortions and the two refocused images are stitched using a simple horizontal translation to match their two common portions. This process can further be totally automatized.

Common refocused images considering all the pixels below the micro-lens are shown at the top of Fig. 11 with and without the prism covering half of the main-lens pupil. The stitched SAI's are visible in Fig. 11c.

Application with a dual-pixel camera

The primary goal of dual-pixel sensor is to perform live auto-focus for video shooting. Recently, with the 5D MarkIV, Canon gives access to the raw dual-pixel. A raw dual-pixel still image contains a first image being the sum of the dual-pixels (i.e. the conventional picture) and a second image being the left SAI. Both images are 30.4M P ix. The right SAI is computed by subtracting the second image to the first one.

The raw images are processed with the following steps: 1/ The raw files are read with [START_REF]Libraw library[END_REF], 2/ The SAI are biascorrected thanks to over-scan areas which are then cropped to keep the responsive pixels, 3/ A white screen is observed with the same lens used with the Raytrix camera, the two SAI's are median filtered, normalized by their maxima, and inverted to define a flux-scale correction which is applied to all captured images. Due to the poor quality of the mainlens, only a central portion of 2000 2 pix is used.

The images captured with and without the prism are shown at the top of Fig. 12 The two SAI's obtained with the prism are stitched to form an extended FOV image (see Fig. 12c). Note that on the stitched image, a ghost image is clearly visible. The sub-aperture cross-talk r is measured masking half of the main-lens while observing a white screen. The cross-talk is as high as r ≈ 50%. A central large mask, masking most of the aperture stop, is placed at the entrance of the main-lens to decrease the cross-talk up to r ≈ 18%, but it remains important. Indeed, we have discovered that the dual-pixel sensor designed by Canon does not produce a sharp sampling of the main-lens pupil. This is certainly done because a sharp sampling would ruin the light efficiency (photons not being collected), and would not improve the auto-focus accuracy. This experiment demonstrates the importance of jointly designing sensors and the corresponding computational algorithms.

Discussion and Conclusion

In the future, ideal quad-pixel sensors could be used in two distinct modes: 1/ the four SAI's are used with typical plenoptic algorithms (e.g. refocus [START_REF] Fiss | Refocusing plenoptic images using depth-adaptive splatting[END_REF], depth-estimation [START_REF] Tao | Depth from shading, defocus, and correspondence using light-field angular coherence[END_REF], partial lens aberration correction [START_REF] Hanrahan | Digital correction of lens aberrations in light field photography[END_REF]) without previous micro-lens center estimation which is a cumbersome task, or 2/ as proposed in this paper, combined with a pyramid lens and stitching the different SAI's to double the FOV as well as the spatial resolution. The pyramid lens would be used as a conversion lens. These two options are especially suitable for smartphones which use a fixed focal lens, and presumably soon a quad-pixel sensor (dual-pixels already integrated).

To complete this study, several technical aspects are to be pointed out: 1/ The SAI's collect photons passing by a quarter disk in the case of a quad-pixel sensor. Thus the bokeh is not a rounded shape as for conventional cameras. Also, we effectively trade the amount of light dedicated to image a scene point against an increase in FOV. 2/ The main-lens has to be designed with an aperture stop located at its entrance and 3/ prisms produce strong chromatic aberrations that degrade the image quality of SAI's, even if achromatic prisms would decrease chromatic aberrations. Finally, using an optical simulation software, the two flat interfaces of the prism could be replaced with quarter spherical interfaces adapted to the main-lens to produce fewer aberrations.

All in all, in this paper we have demonstrated that plenoptic sensors offer a new paradigm for image processing thanks to the ideal pupil sampling they offer. We have presented a novel optical design coupling a pyramidal lens and an ideal plenoptic sensor. The main interest is the capacity to extend the FOV of the captured snapshot. Our design is validated with different experiments including simulated synthetic images and images captured with two different prototypes. In particular, we consider a type-2 plenoptic camera and a Canon DSLR with a dual-pixel sensor, being two available plenoptic cameras on the market. To our knowledge, this is the first time that the scientific community in the field presents experiments with a dual-pixel sensor.

Figure 1 :

 1 Figure 1: Schematic view of a type-1 plenoptic camera.

Figure 2 :

 2 Figure 2: Type-1 plenoptic camera with N = 2. Left: symmetric sampling of the micro-images when φ = 2δ/e. Right: asymmetric sampling of the micro-images when φ = 2δ (considering e = 1).

Figure 3 :

 3 Figure 3: Left: ideal portion of the exit pupil being visible by the SAI S 0,0 when P = 2δ. This portion is constant independently of the sensor coordinate (x, y). Right: portion of the exit pupil visible by S 0,0 for a camera with φ = 2δ. The visible portion is asymmetric and depends on the position of the micro-image on the sensor as well as e.

Figure 4 :

 4 Figure 4: Typical main-lens objective for: single lens reflex camera (left), and smartphone (right).

Figure 5 :

 5 Figure 5: Sensor and main-lens mounted with two prisms.

Figure 6 :

 6 Figure 6: Left: One of the four identical SAI's when there is no prism. Right: Different simulated SAI's obtained with a pyramid lens. The red dot on the left and red square on the right of the test chart indicate respectively the optical axis of the main-lens.

Figure 7 :

 7 Figure 7: Main-lens with three biconvex lenses.

Figure 8 :

 8 Figure 8: Four SAI's showing cross-talk. Ideally, only the bottom-left SAI would capture photons when three quarterpupils are masked but other SAI's are contaminated due to lens aberrations.

Figure 9 :

 9 Figure 9: Crop of a SAI with (right) and without (left) the cross-shaped mask. The cross-talk decreases from r = 13% to r = 3% with the mask. The ghosting artifacts on the left are not visible on the right (both images have the same dynamic range).

Figure 10 :

 10 Figure 10: Type-2 plenoptic camera with P = 4δ = φe Top: Scheme in 1D. Red and green lines indicate the ray beams which exit at the border of the pixel of size δ, converging on two points at plane S and diverging on the mainlens. The colored dashed lines delimit the four projected pixels of one micro-lens into the main-lens pupil. Bottom: Pupil sampling in 2D with the portion of the pupil seen by one pixel in white. Fading regions represent visibility overlapping.

  (a) Full pupil refocus without prism. (b) Full pupil refocus with prism. (c) SAI stitching.

Figure 11 :

 11 Figure 11: Raytrix refocus and stitching.

Figure 12 :

 12 Figure 12: Canon 5D MarkIV refocus and stitching.

Table 2 )

 2 . The shape of an interface is characterized by a central symmetric sag function z(r)

	Thickness Radius Glass Diameter	Conic
	1.82133	3.08301 SK16 3.53718	0.17470
	0.07071	-2.6425	air	3.38069 -18.7316
	0.43551	-2.1581	F5	3.38110 -11.9393
	1.51451	4.72854	air	3.24343	5.13491
	3.79133	3.81723 SK16 4.38871 -4.97297
	1.03669	6.23189	air	4.73167	-1.297

Table 2 :

 2 Description of the six interfaces defining the real main lens.