Loïc Denis
email: loic.denis@univ-st-etienne.fr

Eric Thiébaut
email: eric.thiebaut@univ-lyon1.fr

Ferréol Soulez
email: ferreol.soulez@univ-lyon1.fr

Jean-Marie M Becker
email: jean-marie.becker@univ-st-etienne.fr

Rahul Mourya
email: r.mourya@univ-st-etienne.fr

Fast approximations of shift-variant blur

Keywords: Blur, Deconvolution, Inverse Problems, Image Restoration, PSF

Image deblurring is essential in high resolution imaging, e.g., astronomy, microscopy or computational photography. Shift-invariant blur is fully characterized by a single point-spread-function (PSF). Blurring is then modeled by a convolution, leading to efficient algorithms for blur simulation and removal that rely on fast Fourier transforms. However, in many different contexts, blur cannot be considered constant throughout the field-of-view, and thus necessitates to model variations of the PSF with the location. These models must achieve a trade-off between the accuracy that can be reached with their flexibility, and their computational efficiency.

Several fast approximations of blur have been proposed in the literature. We give a unified presentation of these methods in the light of matrix decompositions of the blurring operator. We establish the connection between different computational tricks that can be found in the literature and the physical sense of corresponding approximations in

Introduction

Image deconvolution is widely used to enhance the resolution, signal-to-noise ratio and contrast of blurred images. This well studied problem is classically stated as the minimization of a suitable cost function composed of a data fidelity and regularization terms [START_REF] Titterington | General structure of regularization procedures in image reconstruction[END_REF]; [START_REF] Demoment | Image reconstruction and restoration: Overview of common estimation structures and problems[END_REF]. The best restored image f reaches a compromise between fidelity to the data g and some priors:

f = arg min f Ψ data (f , g) + µ Ψ prior (f ) , (1)
where Ψ data is a likelihood term derived from the noise statistics and the image formation model, Ψ prior (f ) is a regularization term which has to be minimized to enforce the priors and µ > 0 is a tuning parameter used to set the relative level of the priors.

As far as a small field-of-view is considered, the point spread function (PSF) can be considered as shift-invariant; as a result blurring amounts to a convolution. In many cases, however, blur is space-variant and thus cannot be modeled by a single convolution.

Blur variations across an image can be due to several causes: relative motion between the camera and the scene; moving objects with respect to the background; variable defocusing of non planar scenes with some objects located in front of, or behind the in-focus plane; optical aberrations such as space-variant distortions, vignetting or phase aberrations. Occlusion phenomena occur, e.g., when some objects are in relative motion with respect to the scene's background, or due to parallax when a camera moves during acquisition, or with object masking in scenes displaying varying defocus blurring. There, however, exists a fundamental distinction between configurations with and without occlusion: occlusions arise from the 3-D geometry of the scene, thus a blurred image can not be related to a single planar crisp image. Occlusions are generally handled by compositing layers using their associated alpha matte (i.e., opacity image), see [START_REF] Porter | Compositing digital images[END_REF]. The matting process, i.e., the extraction of foreground and alpha images is well studied for crisp images [START_REF] Chuang | A bayesian approach to digital matting[END_REF]; [START_REF] Wang | Image and video matting: a survey[END_REF]. Deblurring images with space-variant blur involving occlusions requires a simultaneous recovering of the layers of crisp images and transparency channels, which is much more challenging, see [START_REF] Levin | Blind motion deblurring using image statistics[END_REF]; [START_REF] Almeida | Blind deblurring of foregroundbackground images[END_REF]; [START_REF] Chakrabarti | Analyzing spatiallyvarying blur[END_REF].

The focus of this paper is shift-variant blur without occlusion. Several applications with practical interest indeed suffer from smooth blur variations, i.e., occlusionfree phenomena. They include 3-D microscopy, as described in [START_REF] Preza | Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy[END_REF], astronomy with adaptiveoptics correction [START_REF] Cresci | Accounting for the anisoplanatic point spread function in deep wide-field adaptive optics images[END_REF], or wide-field imaging, to name a few. Our focus on occlusion-free phenomena is justified by a methodological concern: PSF variations within the observed field are handled in a completely different way when multi-layer compositing is used compared to singlelayer space-variant blurring. In the case of multiple layers, PSF are defined independently for each layer, and opacities (alpha values) locally govern which layer is visible, and hence, the local PSF. Occlusion-free blur requires modeling smooth PSF variations across the field. Any progress in this direction will be beneficial to matting-based occlusion modeling by handling blur variations within each layer.

While accurate modeling of PSF variations is essential for restoration, iterative deblurring techniques demand fast blurring computation. Blur models must address the two competing objectives of being fast and accurate.

State of the art: When PSFs have very small supports, computation of the blurring operator can be done in the spatial domain, see for example [START_REF] Sorel | Space-variant restoration of images degraded by camera motion blur[END_REF]; [START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF]. PSF with wider supports require blur approximations. The most straightforward approach to shift-variant blur modeling is by decomposition of the observed field into small-enough regions so that the PSF may be considered invariant within each region. Thus these regions can be handled as usual with shift-invariant blur by computing (FFT-based) convolutions, see Sec. 2.2. Such an approach however generates important artifacts at the region boundaries due to PSF discontinuities. In order to reduce these artifacts, Nagy and O'Leary (1998) propose to smooth out the transitions by interpolating between blurred images obtained by convolutions with different PSF. This modification is a strong improvement on piecewise constant PSF modeling for a reasonable additional computational cost (about 4 times more in 2D). This approach lacks physical basis in that it is not related to a natural approximation of PSFs; we show in Sec. 3 that a re-ordering of the operations (namely, weighting before convolving) improves the approximations without increasing the complexity. This formulation has been independently suggested in [START_REF] Gilad | A fast algorithm for convolution integrals with space and time variant kernels[END_REF]; [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF]. We have shown in [START_REF] Denis | Fast model of space-variant blurring and its application to deconvolution in astronomy[END_REF] that it is a natural consequence of PSF interpolation. We recall this result in section Sec. 2.3. A different perspective is followed by [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] and [START_REF] Miraut | Efficient shift-variant image restoration using deformable filtering (part i)[END_REF], see Sec. 2.4, where PSF are decomposed on a few modes computed by principal component analysis. While this paper was under review, we became aware of recent works that use wavelet transforms to efficiently encode the shift-variant blurring operator: [START_REF] Wei | Fast space-varying convolution using matrix source coding with applications to camera stray light reduction[END_REF]; [START_REF] Escande | Numerical computation of spatially varying blur operators a review of existing approaches with a new one[END_REF]. These approaches, not reviewed in details here, offer an appealing trade-off between accuracy and speed, at the cost of a rather intricate relationship between the PSFs and the approximated blurring operator. Finally, rather than modeling shift-variant blur, [START_REF] Maalouf | Fluorescence microscopy three-dimensional depth variant point spread function interpolation using Zernike moments[END_REF] computes several shift-invariant deconvolutions and interpolates the results. Although blur is incorrectly modeled, satisfactory results are shown with this crude approach on some microscopy images.

Our contributions: This paper extends our conference paper [START_REF] Denis | Fast model of space-variant blurring and its application to deconvolution in astronomy[END_REF]. We address the problem of defining fast and accurate models for space variant blur. In a discrete setting, we show that this problem can be recast as a matrix decomposition and approximation problem, thereby providing a unified framework for the description of existing models. We bridge the gap between approaches based on PSF interpolation and decomposition onto PSF modes by defining an optimal local PSF approximation. Both theoretical and practical performance of each model are compared. The formulations of [START_REF] Gilad | A fast algorithm for convolution integrals with space and time variant kernels[END_REF]; [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] appear better grounded and more effective in its approximation of PSF variations than the original method of Nagy and O'Leary (1998) that is still dominant in the literature. The proposed optimal local PSF approximation is shown to improve significantly (typically by an order of magnitude) the approximation error without increasing the computational cost.

Organization of the paper: Section 2 introduces the modeling of shift-variant blur. Existing models cited in the previ-ous paragraphs are described in the light of matrix decompositions. A new approximation with established optimal trade-off between accuracy and computational efficiency is introduced in Sec. 2.5. All these models are compared in section 3 in terms of equivalent PSFs, preservation of physical properties, computational complexity, quality of PSF approximation and image restoration.

Approximations of shift-variant blur

Distortions caused by atmosphere turbulence, object / camera relative motion, or the instrument (limited aperture, optical aberrations) transform the original (crisp) image f into a blurry one, g. A fairly general modeling of this transform takes the form of a Fredholm integral equation of the first kind:

g(r) = h(r, s) f (s) ds , (2)
where h denotes the PSF. The PSF h may be considered as the conditional density p(r|s) describing the probability that a photon entering the system at location s leaves it at location r, see [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF]. In some cases, the PSF is shift-invariant (∀t, h(r, s) = h(r + t, s + t)), i.e., depends only on the difference rs. The right hand side of equation (2) is then a convolution and the system is referred to as isoplanatic.

More generally, the PSF may vary smoothly with the input location s. In order to distinguish true PSF variations from simple shifts of the PSF h(r, s) due to changes in the input location s, it will prove useful in the following to consider un-shifted PSF defined by: k(r, s) = h(r +s, s). Equation (2) can then be rewritten under the form:

g(r) = k(r -s, s) f (s) ds .
(3)

In the general case, the simulation of a blurred image by evaluation of the integral (2) is computationally expensive. As noted by [START_REF] Gilad | A fast algorithm for convolution integrals with space and time variant kernels[END_REF], this evaluation is highly simplified if a separable linear approximation of the kernel is used:

k(r, s) ≈ ∑ p m p (r) w p (s) , (4)
where m p are components of the PSF model and w p are weights depending on the location s. With constant weights w p (s) = w p , the corresponding kernel would be shiftinvariant. By letting the weight w p (s) of each model m p vary with the location s, a shift-variant model is obtained. With this specific approximation, equation (2) simplifies to a sum of convolutions:

g(r) ≈ ∑ p m p (r -s) [w p (s) f (s)] ds ≡ ∑ p m p * (w p f ) (r) ( 5
)
where * is the classical notation for convolution. Equation ( 5) approximates the shift-variant operator as a sum of convolutions of weighted versions of the input image f . The existence of fast algorithms for discrete convolution makes this decomposition extremely useful. It will be central to all approximations that will be studied henceforth.

We now introduce a discrete form for blurring operations modeled by equation ( 2). Approximation of the discrete version of the blurring operator can then be considered from the point of view of matrix decomposition / approximation problems. A linear approximation of the integral in equation ( 2) leads to the following discrete linear model for blur:

g = H • f , (6)
with g ∈ R N the N-pixels blurry image, f ∈ R M the M-pixels crisp image, H ∈ R N×M the blurring operator, and • the matrix product. The matrix H defining the discrete operator is obtained by sampling the continuous operator h at locations (r i ) i=1,..,N and (s j ) j=1,..,M :

H i,j 1≤i≤N 1≤ j≤M = h(r i , s j ) ∆ j , (7)
with ∆ j the elementary volume measure ensuring normalization of H and possible non-uniform sampling of the input field (s j ) j=1,..,M . The j th column H •, j corresponds to the sampled PSF for a point-source located at s j .

Discretization (7) has some limitations. A more accurate discretization of integral (2) can be obtained using generalized sampling theory, as recently described by [START_REF] Chacko | Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography[END_REF]. Let f int be the continuous image defined by using the sequence of discrete coefficients f as weights of a set of basis functions:

f int (s) = ∑ j ϑ int j (s)f j , (8)
with ϑ int j a shifted copy of a certain "mother" basis function ϑ int (e.g., B-splines). Coefficients f are typically chosen so as to minimize the approximation error, i.e., the continuous image f int corresponds to the orthogonal projection of f onto the subspace spanned by basis functions ϑ int j . Moreover, digitization of the blurred image by the sensor involves integration on the sensitive area of the pixel that will be modeled in the following way:

g i = ϑ pix i (r)g(r) dr , (9)
with ϑ pix i a shifted copy of the pixel spatial sensitivity (e.g., indicator function of the sensitive area). Collecting together equations ( 2), ( 8) and ( 9), we get the following set of equations for the f j coefficients:

g i ≈ ϑ pix i (r) h(r, s) ∑ j ϑ int j (s)f j ds dr (10)
Thus the discrete operator H can be defined in a more general form than in equation ( 7):

H i,j 1≤i≤N 1≤ j≤M = ϑ pix i (r)h(r, s) ϑ int j (s) ds dr . (11)
In order to separate changes that modify the PSF from shifts due to the displacement of the point source, we introduce a matrix K that collects all un-shifted PSFs 1 , similarly to the un-shifted kernel k introduced in the continuous case:

K i,j 1≤i≤N 1≤ j≤M = ϑ pix i (r)k(r -s, s) ϑ int j (s) ds dr . ( 12
)
Using separable linear approximations (see equation ( 4))

K i,j ≈ ∑ p m p (i) w p ( j) ↔ K ≈ ∑ p m p • w p (13)
provides decompositions of the shift-variant blurring operator as a sum of convolutions with prior weightings

H ≈ ∑ p conv(m p ) • diag(w p ) , (14)
where conv(m p ) denotes the discrete convolution matrix with kernel m p , diag(w p ) is a diagonal matrix whose diagonal is given by vector w p , and symbol denotes matrix transpose.

Note that if the same kind of approximation were applied to the adjoint operator H , operator H would then be written as a sum of convolutions followed by spatial weighting (rather than preceded by spatial weighting). It would then resemble deformable filtering (a.k.a. steerable filtering), as developed in the seminal works of [START_REF] Freeman | The design and use of steerable filters[END_REF], [START_REF] Greenspan | Overcomplete steerable pyramid filters and rotation invariance[END_REF] and [START_REF] Perona | Deformable kernels for early vision[END_REF]. The key principle of steerable filtering is the ability to compute quickly the (shift-invariant) filtering result for any kernel orientation and size by simple linear combination of few pre-computed (shift-invariant) filtering outputs with well-chosen convolution kernels. The extension of such approaches to shift-variant filtering requires to approximate the rows of discrete operator H (i.e., the integration kernels) rather than its columns (i.e., the PSFs), which may be challenging since only PSF measurements are straightforward. We refer the interested reader to [START_REF] Miraut | Efficient shift-variant image restoration using deformable filtering (part i)[END_REF] 1 this matrix can be considered as a dictionary of all PSFs.

for an extended discussion on shift-variant filtering in the light of deformable filtering2 .

We give in Table 1 an illustration of some K and H matrices on very small-scale 1-D problems with Gaussian PSF. Such representations show how the matrices are structured and how their structure can be exploited to derive powerful decompositions. In contrast to [START_REF] Fish | Scanning singular-valuedecomposition method for restoration of images with spacevariant blur[END_REF], matrices K and H will neither be stored nor explicitly inverted in the following.

Shift-invariant PSF

For a shift-invariant PSF, K is a rank-one matrix with identical columns equal to a single PSF k:

K = k • 1 , (15)
where 1 denotes a vector of convenient length composed of 1. Operator H then corresponds to a discrete convolution.

While circular discrete convolutions are mapped to a simple product by the discrete Fourier transform, discrete convolutions require adequate zero-padding and cropping steps:

H ≡ conv(k) = R •F F F -1 • diag( k) •F F F circular convolution •E , (16)
where E is an expansion operator that adds zeros to the edges of the input signal, R is a restriction operator that truncates the input signal to the size of the visible field of the output image, F F F and F F F -1 are the direct and inverse discrete Fourier transforms, and k is the discrete Fourier transform of the PSF:

k = F F F • E • k . ( 17
)
Applications of operator H to an image f can then be very efficiently computed using fast Fourier transforms (FFT). Correct implementation of discrete convolutions (i.e., without periodization artifacts) is essential in all following shift-variant blur approximations because they rely heavily on discrete convolutions. While circular convolutions applied to shift-invariant blurring generate artifacts located mainly at the image borders, they would yield strong artifacts all over the image when applied to shift-variant models of the form of equation ( 14).

In 1-D problems such as illustrated in the first row of Table 1, convolution matrix H has a Toeplitz structure. For 2-D problems, matrix H is Toeplitz-block Toeplitz and Eq. ( 16) and Eq. ( 17) involve 2-D discrete Fourier transforms. Matrix K is useful to exhibit PSF invariance compared to the more intricate structure of matrix H.

            K =             k • 1 rank one             H = conv(k) discrete convolution piecewise constant (Sec. 2.2)             K =             k 1 • ι 1 +                         conv(k 4 ) •             diag(ϕ 4 ) low rank (Sec. 2.4)             K ≈             u 1 • σ 1 v 1 +             u 2 • σ 2 v 2 optimal global approximation             H ≈               conv(u 1 ) •             diag(σ 1 v 1 ) +               conv(u 2 ) •               diag(σ 2 v 2 ) smooth (Sec. 2.5)             K ≈             c 1 • w 1 +             c 2 • w 2 + • • • +             c 4 • w 4 localized weights optimal local approximation             H ≈             conv(c 1 ) •             diag(w 1 ) + • • • +             conv(c 4 ) •             diag(w 4 ) 1

Piecewise constant PSFs

Let us now turn to the case of PSFs that are shift-invariant within regions defining a partition of the input field. Let P be the number of these regions. There are then P different PSF {k 1 , . . . , k P } and K is now a rank-P matrix:

K = P ∑ p=1 k p • ι p . ( 18
)
As illustrated in Table 1, the j th column of matrix K is equal to PSF k p as long as s j is located in the p th region R p of the input field. The binary weights ι p define the indicator vector:

ι p ( j) = 1 if s j ∈ R p , 0 if s j / ∈ R p . ( 19
)
Within each region, operator H acts like a discrete convolution. It is naturally expandable as a sum of discrete convolution operations with prior weighting in order to restrict each convolution to its specific region:

H = P ∑ p=1 conv(k p ) • diag(ι p ) . (20)
A given column h j of H (i.e., the system response to an impulse at s j ) is thus equal to the corresponding column in the convolution matrix conv(k p ), with p the index of the region containing point s j . The order in which convolution and weighting are applied is essential: weighting first, then convolution. The converse order is used in Nagy and O'Leary (1998) and in many other subsequent works: [START_REF] Calvetti | Restoration of images with spatially variant blur by the gmres method[END_REF]; [START_REF] Nagy | Iterative methods for image deblurring: A Matlab object-oriented approach[END_REF]; [START_REF] Preza | Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy[END_REF]; [START_REF] Bardsley | A computational method for the restoration of images with an unknown, spatiallyvarying blur[END_REF]; [START_REF] Ng | Wavelet restoration of medical pulse-echo ultrasound images in an em framework[END_REF]; [START_REF] Rogers | Strong gravitational lens modeling with spatially variant point-spread functions[END_REF]. The resulting operator is however not equivalent to H, as illustrated in Table 2 and further discussed in section 3. We will denote these approximations under the term image interpolation to emphasize the difference with the PSF interpolation formulation described next.

Smoothly varying PSFs and their local approximation

In the applications discussed in the introduction, PSFs vary smoothly across the field. In such cases, a PSF (e.g., column k j from K) can be well approximated by other PSFs taken at nearby locations (columns k p such that s ps j is small). Let us select P of these columns {k p | p ∈ G P }. Each column of K is then approximated by the weighted sum of P columns out of M (typically with P M):

K ≈ ∑ p∈G P k p • ϕ p , (21)
Interpolation weights ϕ p ( j) are no longer constrained to take binary values. Weights are spatially localized: they are non-zero only on a spatial neighborhood surrounding location s p . The extent of that neighborhood depends on the interpolation order, it corresponds to a square twice the grid step along each dimension for first order (linear) interpolation.

Using this approximation, Eq. ( 14) becomes:

H ≈ ∑ p∈G P conv(k p ) • diag(ϕ p ) . ( 22
)
Weights localization makes decomposition of Eq. ( 22) very suitable for computations: in this way, full-field convolution computations are not necessary since the preceding weighting operation introduces zeros everywhere except on patches with size twice the grid step. This consideration led [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] to suggest to use Eq. ( 22) for shift-variant blur.

The connection with PSF interpolation is however, to the best of our knowledge, not clearly stated in the literature.

We further discuss the consequences of this connection in Sec. 3 where formulation ( 22) is compared with the nonequivalent image interpolation formulation introduced by Nagy and O'Leary:

H ≈ ∑ p∈G P diag(ϕ p ) • conv(k p ) , (23)
in words, convolve then weight.

Low-rank approximation on PSF modes

It is often adequate to consider that PSF variations are well captured by a few number of modes; or in an equivalent way, that K, being a concatenation of PSFs, has low-rank.

A rank-P approximation of matrix K can be expanded as a sum of P rank-one matrices:

K ≈ P ∑ p=1 c p • w p . (24)
The approximation error on K can be expressed as the sum of squared differences between corresponding elements in matrix K and its approximation:

E 2 (A, B) ≡ A -B 2 F = ∑ i, j (A i, j -B i, j ) 2 . ( 25
)
where F designates what is known as Frobenius norm.

In terms of PSF, it is equivalent to the sum of squared Euclidean distances between corresponding PSF (i.e., columns) in K and its approximation. By Eckart-Young theorem, the closest rank-P approximation (with minimum error E 2 ) can be obtained from the singular value decomposition (SVD) of matrix K [START_REF] Martin | The extraordinary svd[END_REF]):

K = U • Σ • V = N ∑ p=1 u p • σ p v p , (26)
Table 2 Comparison of two non-equivalent formulations for shift-variant blur approximation: Approximation is far better using a decomposition that first applies a weighting, then convolves, rather than the converse as proposed originally in Nagy and O'Leary (1998). The curves represent respectively the exact and approximate PSFs with solid lines (resp. dashed lines), sampled at locations s 1 to s 10 . Larger errors are observed in the left panel.

Convolve, then weight Weight, then convolve

Decomposition H ≈ P ∑ p=1             diag(w p ) •             conv(c p ) H ≈ P ∑ p=1             conv(c p ) •             diag(w p ) Recentered PSF approximations k i PSF error (× 4)

Properties interpolates blurred images interpolates PSF preserves positivity preserves positivity -breaks PSF symmetry + preserves PSF symmetry -breaks PSF normalization

+ preserves PSF normalization same computational complexity by retaining only the first P left and right singular vectors, weighted by the corresponding (largest) singular values:

K ≈ P ∑ p=1 u p • σ p v p , (27)
with u p and v p the p-th left and right singular vectors, and σ p the corresponding singular value. In contrast to binary weights of piecewise constant PSF, or localized weights used when interpolating a PSF, components of vectors v p take arbitrary values (positive or negative) and are defined over the whole input field. Vectors u p can no longer be interpreted as PSF (no natural normalization 3 nor positivity), but rather as PSF modes as obtained by principal component analysis, see Table 1.

3 PSF normalization is generally understood in the context of probability density functions: summation of the PSF must be equal to 1 for energy preservation, while it is the (Euclidean) norm of the left singular vectors u p that is taken to be 1.

Decomposition of operator H follows directly from approximation K by reasoning as for Eq. ( 20) and Eq. ( 22):

H ≈ P ∑ p=1 conv(u p ) • diag(σ p v p ) . (28)
Since weights σ p v p are not localized, P full-field convolutions must be computed in this approximation, potentially a large computational load when P 1. The decomposition given by Eq. ( 28) has been proposed in [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] and [START_REF] Miraut | Efficient shift-variant image restoration using deformable filtering (part i)[END_REF].

Proposed optimal local approximation

Low-rank decomposition of Eq. ( 27) is appealing because it is optimal with respect to error E 2 . The corresponding weights are however not localized, increasing the computational cost in proportion with the number of added PSF modes. PSF interpolation approach is preferable in this regard since weights localization prevents the computation of full-field convolutions, saving computational effort especially for small PSF supports. Let us describe an intermediate solution that, while keeping localized weights, is optimal with respect to error E 2 . We define an optimal local approximation of matrix K:

K ≈ P ∑ p=1 c p • w p , (29)
where PSF {c p } P p=1 and weights {w p } P p=1 are optimal solutions in the following minimization problem:

c p , w p P p=1 = arg min {c p ,w p } P p=1 E 2 K, P ∑ p=1 c p • w p . ( 30
)
Without constraint on the support of weights w p , the solution to Eq. ( 30) is given by the truncated SVD (Eq. ( 27)). We however restrict minimization to the support of interpolation weights supp(ϕ p ):

∀p, supp(w p ) ⊂ supp(ϕ p ) , (31)
for a fixed PSF interpolation scheme {ϕ 1 , . . . ,

ϕ P }, see Fig- ure 1.
Eq. ( 30) is a biconvex minimization, as defined by [START_REF] Gorski | Biconvex sets and optimization with biconvex functions: a survey and extensions[END_REF], for which a partial optimum can be found by alternate convex search. Let us describe how each minimization sub-problem is solved.

Computation of optimal PSF: We consider fixed weights w p , gathered into a single matrix W :

W =   | | | w 1 w 2 • • • w P | | |   . ( 32
)
Let C be the matrix whose columns are vectors c p :

C =   | | | c 1 c 2 • • • c P | | |   . ( 33
)
Decompositions ( 29) and ( 30) can be rewritten as a single matrix product:

P ∑ p=1 c p • w p = C • W . (34)
The matrix C of optimal vectors c p , given the weights, is:

C = arg min C K -C • W 2 F = K • W • (W • W ) -1 . ( 35
)
Proof Frobenius norm can be rewritten using the trace function:

K -C • W 2 F = tr[(K -C • W ) • (K -C • W ) ]
whose expansion gives a constant term plus two terms depending upon C, whose derivation with respect to matrix C is:

2 ∂ ∂ C tr(K • W • C ) = 2 K • W , and
∂ ∂ C tr(C • W • W • C ) = 2 C • W • W ,
giving the classical normal equations:

C • W • W = K • W . proving (35).
In most cases a PSF has a support which is much smaller than the observed field. Let L be the number of pixels in the PSF support, N and M the number of pixels of the blurry and original images respectively. Matrix K can be reduced to a L × M matrix by discarding the N -L lines of zeros. Matrix C is then L × P and matrix W is M × P. Eq. ( 35) involves the inversion of a small P × P matrix and computing products involving matrix W . These products are efficiently computed by skipping the many zeros in W due to the restricted supports of the weights (constraint of Eq. ( 31)). [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] approximate local PSF on few PSF modes, the image is thus weighted according to the importance of each mode in the decomposition before convolving with PSF modes; (c) we show in this paper that interpolating PSF leads to the model proposed by [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF], image blocks are first weighted according to the interpolation kernel, then convolved by the PSF; (d) the proposed optimal local approximation of PSF follows the same procedure, with weights and PSF that are chosen so as to minimize approximation error.

Computation of optimal weights: Weights updating requires a minimization of the error E 2 subject to the support constraint of Eq. ( 31). Direct computation of the optimal weights matrix W as done for C in Eq. ( 35) would give non-zero values outside the supports. In fact, only few nonzero weights are used to form the approximation of a given PSF K •, j . These weights are computed independently from weights at other locations. They are obtained by projection of PSF K •, j on the subspace spanned by active PSFs in C, as detailed below.

With the local approximation of Eq. ( 29), element K i, j is approximated by:

K i, j ≈ P ∑ p=1 c p (i) w p ( j) . (36)
A given weight w p ( j) has a global influence on the j th column of K. For any fixed j, only a few w p ( j) are non-zero due to weights localization, typically 2 d in d dimensions for supports corresponding to linear interpolation. Let C j be the sub-matrix obtained by discarding all columns c p of C with index p outside the support of the interpolation weights (i.e., such that ϕ p ( j) = 0). Let the weights inside the support be denoted ω j (i.e., weights allowed to be non-zero), then column K •, j is approximated by C j •ω j . In contrast to Eq. ( 34), the product C j •ω j involves only weights inside the support.

Optimal weights estimation boils down to a least squares formulation:

ω j = arg min ω j K •, j -C j • ω j 2 2 = (C j • C j ) -1 • C j • K •, j . (37)
Since matrix C j is constant within a cell of the interpolation grid, Eq. ( 37) gives the weights of a cell (i.e., intersection of the supports of interpolation kernels; see an illustration of a support drawn with a dashed line in Fig. 1).

C j • C j is a very small 2 d × 2 d matrix (dimension d = 2 for 2D images).
To summarize, minimization of the approximation error is performed by an alternate update of vectors c p with equation (35) and weights w p with equation (37). Starting from PSFs and weights given by the model of linear PSF interpolation described 2.3, each step of the alternate optimization further improves the accuracy of the model until a partial optimum is reached. In practice, only a few iterations are necessary to bring significant improvement (10 iterations where used to get the results illustrated in this paper).

The shift-variant blurring operator H is approximated following the decomposition of Eq. ( 28):

H ≈ P ∑ p=1 conv(c p ) • diag(w p ) . ( 38
)
Optimal vectors c p and optimal weights w p can be computed beforehand (i.e., for a given PSF model H). The complexity of approximation ( 38) is the same as Eq. ( 23).

Figure 2 illustrates how the different models for a smoothly varying blur are applied on an image. Approaches (a), (c) and (d) process small blocks, while the method of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] that approximates PSFs on few modes operates on the whole image. When PSF variations can not be captured accurately using only the first PSF modes, this latter approach is time consuming; approaches based on local processing of image blocks are preferable. The method of image interpolation depicted in Fig. 2(a) is based on the fact that an image with varying blur can be obtained by interpolating between several blurred images computed by convolution with different PSFs. Approximating the blur with PSFs that are (linearly) interpolated from PSFs that are sampled on a grid, as described in section 2.3, leads to the method illustrated in Fig. 2(c). We showed that it is equivalent to the method of [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF]. The proposed optimal local approximation of blur is applied in the same way as the model of Hirsch, using only optimal weights and PSF instead of interpolation weights and sampled PSF.

Comparison of blur approximations

Let us now compare the different blur approximations that have been introduced on the standpoint of preservation of PSF properties (e.g., positivity or symmetry), computational efficiency, approximation error, or deblurring quality.

Equivalent PSF

In the literature devoted to shift-variant blur approximations, the main subject is computational efficiency; most often, the equivalent PSF is not considered. It is however essential to relate a given approximation method with the corresponding approximation in terms of PSF. Some approximations preserve important properties of the original PSF such as normalization, symmetry and positivity. We begin by giving the expression of an equivalent PSF for each model. We then review in section 3.2 the properties preserved by some of the models.

Equivalent PSF for each method: The (centered) PSF k j for a point source located at s j is approximated by an equivalent PSF k j that depends on the model:

k(PSFInterp) j = ∑ p∈G P ϕ p ( j) • k p , (41)
where ϕ p ( j) are interpolation weights;

with the model of image interpolation:

k(ImageInterp) j = ∑ p∈G P diag ϕ → j p • k p , (42)
where ϕ → j p (i) is the interpolation weight at location r i + s j , see derivation in Appendix A; with a model based on decompositions on modes (Sec. 2.4):

k(Modes) j = ∑ p σ p v p ( j) • u p ; ( 43
)
and finally, with the optimal local approximation proposed in Sec. 2.5:

k(OptLoc) j = ∑ p w p ( j) • c p . ( 44
)

Properties

Depending on the approximation method, the equivalent PSF may fulfill some desirable properties that we list now.

Positivity: PSF of systems based on intensity measurements are necessarily positive-valued: ∀(i, j), K i, j ≥ 0.

Equivalent PSF with shift-invariant or piecewise constant models are naturally positive. Models based on PSF interpolation or image interpolation also preserve positivity as long as interpolation weights ϕ are positive. This is in particular the case with linear interpolation. PSF modes are positive-and negative-valued; thus approximations using PSF modes are not guaranteed to preserve positivity. Positivity is not enforced by the proposed local approximation; this can be done by restating criterion (30) under the form of constrained optimization.

Normalization:

The incoming flux at location s j is spread over the support of the PSF h(•, s j ) during the image formation process. If the system is flux-preserving, the sum of PSF coefficients should be equal to one: 1 h j = 1. In practice, there are some losses in the system; more flux may be lost on the border of the field-of-view than in the center, a phenomenon called "vignetting". Correction of this spacevariant attenuation requires an accurate modeling of the total collected flux for each input location, i.e., for all j the sums 1 h j , or equivalently, the sums 1 k j .

A shift-invariant PSF model can only account for a system without space-variant attenuation, as will be illustrated in section 3.5. A piecewise constant PSF model will account for piecewise constant changes of the attenuation; it requires a partitioning into very small regions to model smooth changes. Models based on the interpolation of PSF using positive interpolation weights lead to equivalent PSF with sum 1 k(PSFInterp) j = ∑ p ϕ p ( j) 1 k p . These models then account for variable attenuation by interpolating the attenuation of each sampled PSF k p . In a system with flux conservation, the equivalent PSF conserves the norm of sampled PSF. By contrast, the model based on image interpolation does not preserve summation, even in the case of positive weights where 1 k(ImageInterp)

j = ∑ p,i ϕ → j p (i) • k p (i).
Table 2 summarizes the differences between the widespread model with image interpolation and approaches such as PSF interpolation. It can be observed in the simple case of a Gaussian PSF with standard deviation proportional to location s that the equivalent PSF over-estimates the actual PSF close to s 1 , and under-estimates it near s 10 (so that k(ImageInterp) 1 1 > k 1 1 and k(ImageInterp) 10 1 < k 10 1 ).

Invariance: In many models, the PSFs exhibit some invariance property. Let K be an un-shifted PSF invariant with respect to a linear transform T that preserves un-shifted PSFs: ∀ j, T • k j = k j . Operator T may represent for example a symmetry, or a rotation. Needless to say, shift-invariant and piecewise constant PSF models that approximate matrix K with columns sampled from K preserve the invariance with respect to transform T . The model based on PSF interpolation also preserves this invariance since column k j is approximated by a linear combination of columns sampled from K. By linearity, T • k(PSFInterp) j is then equal to k(PSFInterp) j , i.e., the equivalent un-shifted PSFs are also invariant with respect to transform T . By contrast, the model based on image interpolation that applies varying weights to each sampled PSF does not preserve the invariance property. This can be again noticed on Table 2 where the reference shift-variant PSFs are Gaussian functions and the approximation arising from the image interpolation approach displays some asymmetry, unlike the interpolation model k(PSFInterp) j . The approximation based on a decomposition on PSF modes does preserve invariance properties. Equivalent PSF k(Modes) j are defined in a subspace of the column space of matrix K. Each k(Modes) j could then be re-written as a linear combination of columns of K, and thus is unchanged after application of transform T . Finally, the optimal local approximation method proposed in section 2.5 also preserves invariance properties. The equivalent un-shifted PSF at location s j is defined as a linear combination of vectors c p . These vectors necessarily lie in the column space of matrix K. Otherwise, the rank of K would increase, in contradiction with the optimality of error E 2 .

Computational complexity

The inversion of shift-variant blur requires many evaluations of the approximate blur model. Each approximation discussed in section 2 requires a different computational effort which is analyzed in this section, regardless of the approximation quality that will be studied in the next section. Table 3 summarizes the properties of each model.

Approximating shift-variant blur by an average shiftinvariant PSF is the crudest but fastest method and will be considered as a reference. Let t denote the time required to perform such a convolution for a given image with M pixels and a PSF with a rectangular support of L pixels. The processing time of a piecewise constant PSF model, i.e., constant inside rectangular regions, is also of order t if L M (so that the overhead required to compute values at the outerborder of the regions is negligible). Models based on PSF interpolation have a complexity that depends on the number of dimensions d along which PSFs vary (typically, d = 1, 2 or 3) and on the interpolation order o. With first-order interpolation and 2D shift-variant blur, PSFs are approximated by bi-linear interpolation; there are 2 2 non-zero terms in the sum of equation ( 22). More generally, there are (o+1) d nonzero terms and if outer-border computations are negligible (the support of weights ϕ p being large compared to the support of the PSF), the total computation time is ≈ t ×(o+1) d .

A more accurate expression for complexity can be given. Let us consider a rectangular grid of P PSF sampled uniformly. The number of pixels of the crisp image is M, the number of pixels in the support of the PSF is L, and the number of pixels inside a grid cell is on average M/P. Due to the overlap of interpolation weights, each convolution in equation ( 22) is carried on (o + 1) d cells. These convolutions can be computed using fast Fourier transforms (FFT) with an appropriate padding with L 1/d zeros along each dimension. FFTs are then computed on areas of about [(o + 1) × (M/P) 1/d + L 1/d ] d pixels. Since the whole operator involves P such computations, the total complexity is of order

P d [(o + 1) × (M/P) 1/d + L 1/d ] d log[(o + 1) × (M/P) 1/d + L 1/d ].
If the PSF support is much smaller than The method based on image interpolation and ours have the same complexity as the method described for PSF interpolation since convolutions are computed on areas that have similar sizes. By contrast, the decomposition on PSF modes described by [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] does not enforce localization of weights. Thus, each of the P convolutions must be computed on the full image support, which is much more costly.

The size of the regions affects the efficiency of FFTs and other operations, which is not visible from the asymptotic computational complexity analysis. Smaller regions better fit into the different levels of memory cache; thus, it may be beneficial to use a finer grid of PSFs, as illustrated in table 4 case (a); larger PSFs involve much larger border effects and rapidly increasing computational times when the grid is refined, as illustrated in case (b).

Implementations of the approximations of shift-variant blur can benefit from parallelism at different levels: SIMD instructions for FFT and multiplications and additions, multi-threaded FFT supported by FFT libraries such as FFTW described in [START_REF] Frigo | A Fast Fourier Transform compiler[END_REF], and coarse-grained parallelism where each term of the sum is computed in parallel in a different thread. In our implementation, we used both SIMD instructions and coarse-grained parallelism. Table 5 shows that our code scales well with the number of CPU cores: an average speedup of 5.3 is measured with 6 cores in our experiments. Using the same parallelization strategy for the image interpolation approach and the optimal method proposed in section 2.5, we reach similar timings.

Approximation error

The computational complexity should not be dissociated from the approximation error achievable by each model.

Let us first discuss about the ability of each model to exactly match the reference shift-variant operator when enough terms are used for the approximation. The piecewise constant PSF model matches the reference operator H when the number of terms P equals the number of input pixels M. Each region is then restricted to a single pixel and indicator vectors ι j are elements e j from the canonical basis (i.e., for all (i, j), e j (i) = 0 when i = j and e j (i) = 1 when i = j). Similarly, using interpolation weights ϕ j restricted to a single pixel with ϕ j = e j gives a PSF interpolation model that exactly matches the operator H. By contrast, the model based on image interpolation with the same interpolation weights produces an approximation of operator H with an error that is bounded from below (with the consequence of a systematic irreducible error). Setting the interpolation weight ϕ j equal to the values e j in the expression of the equivalent PSF given in equation ( 42), one obtains: lim

ϕ j →e j k(ImageInterp) j (i) = k(r i , s q r i + s j )∆ j k q (i) = k(r i , s j )∆ j k j (i) . (45)
In the extreme case of a grid of PSF with the same density as the pixel grid, i.e., with as many terms as pixels in the image for the model based on image interpolation, the equivalent PSF does not correspond to the reference PSF. Only the central value (r i = 0) of the equivalent PSF matches that of the reference PSF. Other values are sampled from PSFs located increasingly farther as the output location r i in the PSF moves away from the PSF center. The model based on image interpolation offers a similar computational cost as PSF interpolation or optimal localized PSF/weights without reaching perfect approximation with regions as small as a single pixel. This is a serious reason for this model to be disregarded.

Finally, the model based on a decomposition on PSF modes provides an exact representation of the original operator H as long as the number P of terms is at least equal to the rank of H (at most min(M, N)). In the case of optimal localized weights described in section 2.5, the approximation error is zero in the limit case already considered. The Table 6 Zernike coefficients used in the "two phase-screens" shiftvariant simulation. Only non-zero coefficients are displayed.

aberration name

Zernike coefficients defocus a 4 = 0.3 a 4 = 0.1 third order astigmatism a 6 = 1.4 a 6 =-1.4 third order spherical a 11 = 0.1 a 11 =-0.02 fifth order coma a 16 = 0.05 a 16 = 0 fifth order coma a 17 = 0.02 a 17 = 0 fifth order spherical a 22 =-0.5 a 22 = 0.5 choice c j = k j and w j = e j for all j corresponding to the PSF interpolation model indeed leads to a zero approximation error (E 2 = 0) and is thus a solution to the optimization problem defined in equation ( 30).

In order to compare the different approximations of shift-variant blur on a realistic case, we modeled an optical system with wavefront aberrations located in two distinct pupil planes, see Figure 3. These aberrations may come from inhomogeneity of optical indices in the turbulent layers of atmosphere crossed by light in ground-based astronomical observations, or from optical aberrations of a multi-lens system. Depending on the incoming direction of the plane wave created by a point source at infinity, the optical aberrations do not sum up the same way, resulting in a shift-variant blur. Partial pupil masking for large incident angles creates an attenuation effect called vignetting. Variations of the incidence angle of the incoming wave lead to changing tip-tilt aberrations resulting in different shifts of the PSF centers. Using Fraunhofer diffraction theory, [START_REF] Goodman | Introduction to Fourier optics[END_REF], the PSF writes:

h(r, s) ∝ Γ (ρ) exp j P(ρ) + 2π r • ρ λ dρ 2 , (46)
with Γ (ρ) the modulus of the complex pupil function:

Γ (ρ) = Π (ρ)Π (ρ -s/s max ) , (47)
and with P(ρ) the aberrated phase:

P(ρ) = 2π λ ref λ ∑ k a k Z k (ρ) + a k Z k (ρ -s/s max ) , (48)
where Z k is the k-th unnormalized Zernike polynomial as defined in [START_REF] Mahajan | Zernike circle polynomials and optical aberrations of systems with circular pupils[END_REF], a k and a k are the Zernike coefficients of the first and second phase screens, s and r are angular coordinates in radians, j stands for the imaginary unit, Π (ρ) is the indicator function of a disk of unit radius centered at the origin: Π (ρ) equals one if ρ 2 ≤ 1, otherwise Π (ρ) is zero, λ is the optical wavelength and λ ref = 546 nm is a reference wavelength. Coefficients {a k } and {a k } used in our numerical experiments are given in table 6. Figure 5 displays approximation error maps for the different shift-variant models in two cases: a coarse approximation using a 5 × 4 grid of PSFs, and a refined approximation with a 20 × 16 grid. For the model of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF], we set the number of modes to 5 and 13 which represent a computational effort similar to the use of PSF grids. Unsurprisingly, using a shift-invariant PSF or piecewise constant PSF lead to the worst errors. PSF interpolation and image interpolation approaches are comparable when using a coarse grid, while the error is twice as much with image interpolation when the grid is refined (see subfigures (i) and (j)). Decomposition on PSF modes using the method of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] is penalized by the increase of the computational cost in proportion to the number of modes (errors with 13 modes are over 4 times larger than with PSF interpolation on a 16 × 20 grid).Such an observation is however application-specific since the approximation error with the method of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] depends on the spectrum of singular values of matrix K. Our simulation yields small shifts of the PSF (i.e., geometrical distortions) that reduce strongly the redundancy of the PSF (i.e., the decrease of singular values is slower than it would be with centered PSF). Finally, the highest accuracy is reached using the proposed optimization of localized weights and PSF. Approximation errors are reduced by an order of magnitude without increasing the complexity required to evaluate the operator.

The gain obtained by using weights and PSF refinement is further illustrated on Figure 6. This figure displays the evolution of the approximation error computed on the PSF or on a blurred image as a function of the density of the grid of PSF (solid lines). Approximation errors decrease when the grid is refined, except in some cases where a slight increase of the grid size gives a less favorable sampling of the PSF. From the normalized approximation errors given in log scale on the left-hand side y-axis one sees that PSF and weights optimization improve the approximation by about an order of magnitude. Figure 6 also displays the evolution of the computation times scaled as multiples of a full-size convolution. Increasing the grid size degrades the performance (except for some grid sizes that offer better FFT algorithms or better fit into memory caches). A trade-off must then be reached between the accuracy of the approximation and the time required to apply the operator. Thanks to the improved accuracy of the model with optimized PSF and weights, a given approximation error can be reached at a much lower computational cost. Fig. 6 Trade-off between reduction of the approximation error and increase of the computational cost.

Deblurring performance

The blur models discussed in this paper are ultimately used within image restoration methods. In order to emphasize the impact of data modeling (i.e., the choice of a specific blur model), we considered numerical simulations in which degraded images were computed using a different PSF at each pixel of the field of view. In the first numerical experiments, no random noise was added, so the discrepancy between the simulated image and the output of a given blur model can only reflect modeling errors. The addition of noise tends to alleviate blur model differences. We used the same physically-grounded PSF model as described in Figure 3 and in the previous paragraph. All our image restorations were carried out by solving the same (convex) MAP optimization problem (1), i.e., only the blur model H and the regularization weight µ change. Thus, differences in restoration quality can be directly ascribed to the blur model used. Assuming i.i.d. Gaussian noise, the data fidelity term becomes a simple least squares criterion:

Ψ data (f ) = H • f -g 2 2 . ( 49
)
In order to induce smoothness while preserving edges, we used an hyperbolic approximation of the total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]) for the regularization of our images:

Ψ prior (f ) = ∑ j ∇ j • f 2 2 + ε 2 , (50)
with ∇ j a finite difference operator which approximates the spatial gradient ∇ f (s j ) around j th position s j and ε > 0 a relaxation parameter to have a criterion that is differentiable whatever f . Thanks to this hyperbolic approximation, the image reconstruction problem (1) is a smooth convex minimization problem that can be solved using standard methods for large-scale optimization. We used a limited memory iterative quasi-Newton method4 [START_REF] Nocedal | Updating quasi-Newton matrices with limited storage[END_REF]). Another advantage of introducing parameter ε is to achieve an 2 -1 behavior and avoid the classical staircase artifacts created by total variation minimization, [START_REF] Mugnier | Myopic deconvolution from wave-front sensing[END_REF]. Choosing the best tuning parameter µ in equation ( 1) is the subject of intensive researches. For our purposes, since we knew the ground truth image, we simply tuned µ so as to maximize the PSNR. Note that many alternative regularization terms could be considered, such as recent models based on patch statistics: [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]; [START_REF] Chen | Insights into analysis operator learning: A view from higher-order filter-based mrf model[END_REF].

To compare the performance of the shift-variant blur models described in this paper, we first consider the restoration of a simple test image made of white stripes on a black background. The simulated blurred version is shown in figure 7(a). Below each image a line profile is plotted corresponding to the red horizontal line. Blur is stronger in the periphery of the image compared to its center. Restoration of this test image using an invariant blur model and the central PSF produces the result shown in figure 7(b). The resolution in the center of the field is slightly improved but most of the blur remains in the periphery. Shift-variant models of Nagy and O'Leary (1998) (figure 7(c)) and bi-linear interpolation of PSFs (figure 7(d)) give comparable performance when a coarse 4 × 5 grid is used. The model of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] gives better performance when 5 modes are selected (figure 7(e)), probably because the error in PSF modeling is far more homogeneous in the field compared to the two previous models (see figure 5(e) compared to 5(cd)). It can be noted that the location of the white stripes corresponds to areas where the errors in PSF modeling are the highest for shift-variant models with a 4 × 5 grid. The proposed blur model with optimal PSF and weights gives a far better restoration, as can be visually assessed on figure 7(f) and in terms of PSNR (9dB improvement compared to the shift-variant model with bi-linear interpolation of weights, with identical computational complexity once the optimal PSFs and weights are computed). Restoration results with 16 × 20 refined grids or more modes are shown on figure 7(g-j). The difference between image interpolation and PSF interpolation is clearly visible with this denser sampling of PSFs (a 3dB improvement is observed in favor of PSF interpolation). The model of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] offers a weaker improvement when the number of selected modes changes from 5 to 13 modes than other models with similar computational costs. The proposed model with optimized PSFs and modes reaches the best restoration performance. It can be observed that the value of the optimal weight µ of the smoothness term is smaller when the modeling errors are smaller, i.e., the model of the data is more reliable. The values of µ reported in figure 7 show a reduction by a factor 150 when replacing the shift-invariant model with the most accurate shift-variant model (i.e., with optimized weights and PSFs on a 16 × 20 grid).

We then consider the restoration of a natural image displaying different textured areas, figure 8(a). Since the image is in color, we considered a different PSF model for each channel. The phase aberration model described in figure 3 is adapted to different wavelengths by scaling. The resulting blurred image is shown in figure 8(b). We considered that the sensor covers only the area within the red box superimposed on all images. Let γ be the indicator of that area (i.e., γ k = 1 iff pixel k is in the area covered by the sensor, otherwise γ k = 0). We express our data-fitting term as:

Ψ data (f ) = (H • f -g) diag(γ) (H • f -g) . (51)
In order to correctly model the image borders, we reconstruct an image that extends outside the actually observed area, see [START_REF] Soulez | Inverse problem approach in particle digital holography: out-of-field particle detection made possible[END_REF]; [START_REF] Reeves | Fast image restoration without boundary artifacts[END_REF]; [START_REF] Matakos | Accelerated edge-preserving image restoration without boundary artifacts[END_REF]. We chose to apply separately the relaxed total variation regularization to each color channel, instead of using the color total variation introduced in Blomgren and Chan (1998). Thus, artifacts introduced by each blur model are better identified (i.e., the regularization term does not compensate for inaccuracies of the blur model). The image is 400 × 320 pixels, PSFs are 51 × 51 pixels, ε is set to 10 (values in each color channel are in the range [0, 255]), and 500 iterations of the quasi-Newton method are performed (it has been checked in practice that 300 iterations were generally sufficient and that 2000 iterations did not bring noticeable improvement). Restoration using a shift-invariant model gives the image shown in Figure 9(a). Several artifacts can be observed: the original flux is not restored; geometrical distortions are not corrected; and the original resolution is recovered only in the center of the field of view. Using a shift-variant blur model such as PSF interpolation leads to restored images with a more homogeneous resolution. When a coarse PSF grid is used, PSF interpolation and image interpolation produce very similar results, thus only the result for PSF interpolation is given in figure 9(b). The decomposition on few PSF modes fails to correctly model PSF tilt (i.e., decentering). The restored image shown in figure 9(c) suffers from strong color artifacts. The decomposition on PSF modes is not the best method for a simultaneous correction of shift-variant blur and geometrical distortion (i.e., variable PSF decentering). The model with optimized PSF and weights on a small 4 × 5 grid leads to the best PSNR among coarse models and a significantly improved resolution (e.g., region B in the zooms is much better recovered). However, abrupt changes of the PSF model at the boundary between two regions of the grid produce annoying artifacts. These artifacts can be reduced either by considering a refined grid, or by increasing the regularization weight at the cost of a resolution loss. They can be considered as a result of a misuse of the refined model, the grid being too coarse with respect to PSF variations. Figure 9(e-h) gives deblurring results with fine blur models of comparable computational complexity. With a denser grid of PSFs, the difference between deblurring with a model based on image interpolation (figure 9(e)) and PSF interpolation (figure 9(f)) is now noticeable in terms of PSNR (+0.25dB). When going from 5 to 13 modes, the results obtained with the model of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] improve notably. Within a comparable time budget, this model however gives less accurate restorations (PSNR is 2dB smaller). The proposed optimization of weights and PSF again leads to the best results in terms of PSNR. Artifacts at the boundaries of regions that could be observed with a coarse grid disappear with finer grids (compare figures 9(d) and 9(h)).

Table 7 provides further comparison of blur models based on the PSNR improvement brought by shift-variant blur models compared to a shift-invariant model based on the central PSF. Five natural images and three levels of noise are considered. Noise levels are given in terms of BSNR, i.e., the ratio between the noiseless blurred image range and the noise standard deviation, in log scale. In addition to Barbara image, we chose four images from Berkeley's segmentation dataset5 . Some of the deblurring results of table 7, denoted by a star, are displayed in figure 10. Like in the previous experiments, we selected the number of modes to get comparable computational times with all blur model. These deblurring results show a significant improvement for all models that model blur variations compared to the shift-invariant model. While the use of an accurate model brings small improvements for moderate SNRs, deblurring is strongly improved by accurate blur modeling in the high SNR regime. Interestingly, improving the SNR brings almost no improvement when using a shift-invariant blur model, for the blur under consideration. This means that in our setting modeling errors dominate stochastic errors from noise, so that only a model refinement can improve the deblurring results. Blur models based on a refined grid or a larger number of PSF modes consistently improve on coarser models. Among each group of models with similar computational complexity, the proposed locally optimal approximation provides the best deblurring performance, sometimes on par with the approach based on PSF interpolation. PSF interpolation and image interpolation be- Table 7 Restoration performance in terms of PSNR gain in dB with respect to shift-invariant deblurring. The noise level used in our simulations is given in BSNR (see text). The images for the numbers followed by a star are displayed in Fig. 10. have similarly when a coarse grid is considered, or for lower values of the SNR. The improved accuracy brought by PSF interpolation models is visible in the high SNR regime and with refined grids (size of the grid should be comparable or smaller than the PSF size). By close inspection of the images in figure 10, block artifacts are visible with the proposed locally optimal method only in the case of very high SNR and coarse grid (sixth row of images, last column).

Conclusions

We have provided a unified presentation of several shiftvariant blur models that co-exist in the literature. Some of these models had earlier been introduced merely as computational tricks to produce a shift-variant blur effect, within a time lapse comparable to few convolutions. Detailed analysis of how those models relate to particular PSF approximations show that some are more sound, physically speak- ing, thus leading to improved performance for shift-variant deblurring. The most natural model follows from assuming that the PSF between two given locations can be interpolated, [START_REF] Denis | Fast model of space-variant blurring and its application to deconvolution in astronomy[END_REF]. It leads to a numerical implementation as efficient as the more widespread model based on image interpolation due to Nagy and O'Leary (1998), but with a better grounded approximation yielding systematic improvement of the deblurring quality. Incidentally, this PSF interpolation model has the same mathematical formulation as in [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF]. These blur models can be refined in order to achieve high accuracy, at the cost of increased computational cost. We proposed to seek, within a given computational budget, the best possible model in terms of approximation quality. We modified the model based on PSF interpolation by considering optimal interpolation weights and PSF samples. This quite simple method improves the RMS modeling error of PSFs by as much as an order of magnitude. When applied to deblurring, this model consistently provides improved restorations, as quantified by PSNR values. Some block artifacts can appear if the grid of PSF is too coarse. These artifacts are due to some discontinuities of PSF approximations between neighboring blocks. They are visible only for high signal to noise ratios and disappear either by increasing the regularization or by refining the grid.

The cost of accurate models can be modest (an increase of computational time typically between a factor 4 to 10 in 2D) compared to convolution by a stationary PSF. The computational complexity should no longer be an obstacle to account for blur variations within the field of view.

Modeling a shift-variant PSF is necessary but not sufficient. The next important task is to properly calibrate this PSF. In practice, this is considerably more difficult than for a shift-invariant PSF because it requires the ability to put reference sources at any spot. The methodology proposed in [START_REF] Delbracio | The nonparametric sub-pixel local point spread function estimation is a well posed problem[END_REF] seems a promising approach to estimate shift-variant PSFs. A more general method would be to follow a blind deconvolution approach which is able to estimate the deblurred image and the blur from a single blurred image, see e.g., [START_REF] Ayers | Iterative blind deconvolution method and its applications[END_REF]; Campisi and Egiazarian (2007); [START_REF] Soulez | Blind deconvolution of 3d data in wide field fluorescence microscopy[END_REF] and [START_REF] Hadj | Blind restoration of confocal microscopy images in presence of a depth-variant blur and poisson noise[END_REF]. In some sense, the PSF and weight optimization strategy that we used in order to improve the quality of the model of the shift-variant blur is already a step in that direction.

A Derivation of the equivalent PSF for the model of Nagy and O'Leary (1998) Nagy and O'Leary (1998) proposed to approximate the shift-variant blurring operator H by interpolating between the result of several shift-invariant blur:

H ≈ ∑ p∈G P diag(ϕ p ) • conv(k p ) . (52)
Compared to expression (22) for PSF interpolation, the order of the weighting and blurring operations is exchanged. The right-hand side of equation ( 52) does not correspond to blurring with an equivalent PSF that is interpolated from a few sampled PSF. The equivalent un-shifted PSF k j at location s j , as defined in section 2, is: k j (i) = k(r i , s j )∆ j = h(r i + s j , s j )∆ j . From equation ( 52), we get:

h(r i , s j ) = ∑ p∈G P ϕ p (r i ) • k p (r i -s j ) . ( 53
)
The equivalent PSF is then:

k j (i) = h(r i + s j , s j )∆ j = ∑ p∈G P ϕ p (r i + s j ) • k p (r i + s j -s j )∆ j = ∑ p∈G P ϕ → j p (i) • k p (i) , (54)
where ϕ → j p designates the interpolation weights ϕ p shifted by the vector s j . Using vector notation, the equivalent PSF can be written as in equation ( 42) that we recall here:

k(ImageInterp) j = ∑ p∈G P diag ϕ → j p • k p .

Fig. 1

 1 Fig. 1 2D illustration of shift-variant blur approximations: one of the convolution kernels and weights used with (a) piecewise constant PSF, Sec. 2.2; (b) interpolated PSF, Sec. 2.3; (c) decomposition on PSF modes, Sec. 2.4; (d) proposed optimal local approximation, Sec. 2.5. In (b) and (d), weights are non-zero only on a small region drawn with a dashed line (for visualization purposes, a 4 × 4 grid of PSF is considered here, weights are much more localized with denser PSF grids).

Fig. 2

 2 Fig.2Shift-variant blur applied to an image with 4 different models: (a) the model ofNagy and O'Leary (1998) first convolves image regions with different PSF and then interpolates the blurry results; (b)[START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] approximate local PSF on few PSF modes, the image is thus weighted according to the importance of each mode in the decomposition before convolving with PSF modes; (c) we show in this paper that interpolating PSF leads to the model proposed by[START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF], image blocks are first weighted according to the interpolation kernel, then convolved by the PSF; (d) the proposed optimal local approximation of PSF follows the same procedure, with weights and PSF that are chosen so as to minimize approximation error.

=

  k p , with p such that ι p ( j) = 1 ; (40) with a model based on PSF interpolation (Sec. 2.3):

Fig. 3

 3 Fig.3Shift-variant blurring due to phase aberrations and vignetting: an incident plane wave is distorted according to phase aberrations of each lens. In the center of the field (b), most phase aberrations are compensated. Farther from the center of the field, part of the incoming wave front is masked (vignetting), and phase aberrations create a tilt and larger spread of the PSF.

Fig. 4

 4 Fig. 4 Grid of PSFs generated from the shift-variant model depicted in Fig.3. Contrast is inverted in order to improve its visualization.

Fig. 5

 5 Fig.5Approximation errors of a shift-variant blur. Error maps shown give the RMS errors between the true PSF and its approximation depending on the location of the PSFs.

  a full-size convolution) normalized rms error size of the grid of PSF used in the fast approximation

Fig. 7

 7 Fig. 7 Restoration of a resolution target degraded by shift-variant blur: (a) degraded image, (b) single-PSF deblurring,(c-f) deblurring with shift-variant PSF models of comparable computational complexity (coarse models), (g-j) deblurring with shift-variant PSF models of comparable computational complexity (fine models). A line profile along the red line indicated by the symbols and is drawn below each image.

Fig. 8

 8 Fig. 8 Image used to evaluate shift-variant deblurring of a color image: (a) original image; (b) simulated blurry image with a shift-variant chromatic PSF.

Fig. 9

 9 Fig. 9 Restoration of a color image degraded by chromatic and shift-variant blur: (a) deconvolution with a shift-invariant model; (b-d) deblurring with shift-variant PSF models of comparable computational complexity (coarse models).

Fig. 9 (

 9 Fig. 9 (cont.) Restoration of a color image degraded by chromatic and shift-variant blur: (e-h) deblurring with shift-variant PSF models of comparable computational complexity (fine models).

  Fig. 10 Illustration of some of the deblurring results from the numerical experiment of Tab. 7.

Table 1

 1 Decomposition schemes of PSF K and blurring operator H in different scenarios of PSF variations. Columns of matrices K and H represent PSFs. Matrix entries are represented as dots whose area is proportional to their magnitude. Negative values are displayed as circles.

	PSF variations	matrix	decomposition	properties
	Note: invariant
	(Sec. 2.1)

Table 3

 3 Summary of the main properties of shift-variant blur models (P is the number of terms in the approximation) is small compared to the size of the regions; for approximations involving the 4 nearest PSFs References: [A] Maalouf et al (2011); [B] Nagy and O'Leary (1998); [C] Hirsch et al (2010); [D][START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF]

	Method	Reference	Assumptions	Properties	Complexity
					(convolutions)
	interpolate deconvolution results	[A]	slow PSF variations	-no shift-variant PSF model	≈ P
	piecewise constant PSF	Sec. 2.2	large isoplanatic regions	-strong boundary artifacts	≈ 1
	convolve, then apply linear weighting [B], Eq. (23)	smooth PSF variations	+ preserves PSF positivity	≈ 4 in 2D
				+ interpolates PSF,
	use linear weighting, then convolve	[C], Sec. 2.3	smooth PSF variations	+ preserves PSF positivity, normalization and symmetry	≈ 4 in 2D
	decompose on PSF modes	[D], Sec. 2.4 PSF captured by few modes + optimal global approximation	P
	use optimal weighting, then convolve	Sec. 2.5	smooth PSF variations	+ optimal local approximation	≈ 4 in 2D
	if PSF support

Table 4

 4 Computation time relative to a simple convolution for an approximation of operator H based on linear interpolation of PSF (Sec. 2.3): (a) 512 2 pixels image, with 31 × 31 pixels PSFs; (b) 1000 2 pixels image, with 101 × 101 pixels PSFs. On our 3.3 GHz Intel Xeon Processor, a simple convolution takes 19 ms for a 512 2 pixels image (case a) and 43 ms 1000 2 pixels image (case b).

	PSF grid: 5 × 5 10 × 10 20 × 20
	case (a):	3.7	2.4	3.2
	case (b):	4.2	15	39
	each cell of the grid of PSF (L M/P), the complexity is of
	order (o + 1) d × M log M, i.e., it corresponds to (o + 1) d full-size convolutions (with first-order interpolation, it yields a
	multiplication by 2 d , where d = 1, 2 or 3 is the dimension).
	With very big PSF supports (L		M/P), the complexity
	rises to P full-size convolutions.

Table 5

 5 Average computation time for an approximation of operator H based on linear interpolation of PSF (Sec. 2.3): (a) 512 2 pixels image, with 31 × 31 pixels PSF; (b) 1000 2 pixels image, with 101 × 101 pixels PSF (Intel Xeon Processor with 6 cores). The speedup factor relative to a single-threaded computation is placed between brackets.

	nb. of cores:	1	2	4	6
	PSF grid	time in ms (speedup gain w.r.t. one core)
	5 × 5 (a) 10 × 10 20 × 20	71 46 60	37 (1.9) 22 (2.1) 30 (2.0)	19 (3.7) 11 (4.2) 15 (4.0)	14 (5.1) 8 (5.8) 11 (5.5)
	182 636 325 (2.0) 170 (3.7) 117 (5.4) 97 (1.9) 57 (3.2) 47 (3.9) 20 × 20 1680 850 (2.0) 435 (3.9) 293 (5.7) 5 × 5 (b) 10 × 10

in particular, the adjoint operator H of operator H defined in equation (14) corresponds to a deformable filter

the source code can be downloaded at: http://cral. univ-lyon1.fr/labo/perso/eric.thiebaut/?Software/ OptimPack

http://www.eecs.berkeley.edu/Research/Projects/ CS/vision/bsds/

This work has been supported by project MiTiV funded by the French National Research Agency (ANR DEFI 09-EMER-008-01). It has been performed in part within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). Rahul Mourya acknowledges a PhD grant funded by the Région Rhône-Alpes.


Clement Poncelet

Florent Jacquemard

Compilation of the Intermediate Representation V1

Keywords: Model-based testing, automatic model construction, Antescofo

The IR has the form of an executable code modeling the expected Antescofo's behavior on a given score. We present here a simplified version of Antescofo's IR suitable for our presentation, leaving features such as thread creation, conditional branching, and variable handling outside of the scope of this paper. In this paper, we define the construction of the first version of the IR from a given Antescofo DSL score.

Compilation de la Représentation Intermédiaire v1

Résumé : La Représentation intermédiaire est sous forme d'un code exécutable modèlisant le comportement attendu d'Antescofo pour une partition donnée. Nours présentons ici une version simplifiée de cette représentation, mettant de côté la création de thread, les branchements conditionnels et la gestion de variable. Dans ce papier, nous définissons la construction de la première version de la représentation intermédiaire depuis la partition Antescofo.

Mots-clés : test basé sur model, construction de modèle automatique, Antescofo

introduction

1.1 The score-based IMS Antescofo representation

The IR construction is based on the parsing of an Antescofo mixed score abstract syntax tree. We represent this AST (the handled part) as described Figure 1.

Let O be a set of output messages (also called action symbols and denoted a) which can be emitted by the system and let I be a set of event symbols (denoted e) to be detected by the LM (i.e. positions in score). An action is a term act(d, s, al ) where d is the delay before starting the action, s is either an atom in O or a finite sequence of actions (such a sequence is called a group), and al is a list of attributes. A mixed score is a finite sequence of input events of the form evt(e, d, s) where e ∈ I, d is a duration and s is the top-level group triggered by e. Sequences are denoted with square brackets [, ] and the empty sequence is [ ].

A sequence of actions is a suite of actions that can be another group or an atomic action. An atomic action can be a message or a kill action. Messages and groups can have attributes for their synchronization strategy and/or error management (@loose @global by default). We distinct an atomic action of a group by the presence or not of an actions sequence in al.

Model construction

The construction is done through a compiler, which visits the AST (design pattern visitor) and constructs the corresponding network of machines. Currently we represent graphically these machines with extended FSMs.

FSM representation

During the parsing, part of machines are built and merged with other (previous and next parts).

For this purpose we use FSMs with two distinguished sequences of locations called providers and seekers. They are marked with i•, resp. i, where i is an index in the sequence of providers, resp. seekers. The operator G 1 +G 2 takes two FSMs with the same number of seekers for G 1 as providers for G 2 , and returns one FSM G obtained by a union of G 1 and G 2 where moreover, every seeker of G 1 is merged with the provider of G 2 with the same index. The sequence of providers (resp. seekers) of G is the sequence of providers of G 1 (resp. sequence of seekers of G 2 ), see the Figure below.

G 1 0 1 2 τ 1 ? τ 1 ? + G 2 1 2 i 1 d τ2! 2 τ 2 ! = G 0 1 1 τ 1 ? τ 1 ? 2 1 d τ2! 2 τ 2!

Inference rules

Parsing the Score Inference rules are used to explain the recursive behavior of each case, during the parsing phases. Keep in brain that a rule presents in its top what you need to do (one or more sub-functions, or itself recursively) for doing a step of the function at its bottom part. A function is represented with its inputs arguments, its type name ( ), the current AST parsed part and the IR returned. Inference rules describe the parsing and model construction of the compiler. It is depicted as:

input AST : IR • • • input AST : F inalIR ex
with F inalIR = IR + rule step and < input > = acc, sc, e, an accumulator of time elapsed since the beginning of processing of the group, a score sc and a current event e (called the next closest event). Generally the parsed AST ( AST ) is an event for rules (EG,...) which contains a group (otherwise nothing is done and the event is simply past), or an actions sequence with a current atomic action or group parsed. In the latter case rules MSG or GRP respectively is executed.

Beginning

We start with two rules, which parse the entire AST ( all ). The first rule expresses the empty case returning an empty FSM called A ∅ . The other starts by creating three parallelized machines (expressed with the operator ||), E for the environment, S for the Specification composing by P the proxy and G the sub-network of machines modeling Antescofo score-specific behaviors. We let M the entire IR-model composing by the environment E and the specification S and we define the operator || as a current parallelization of two states machines.

Environment and proxy

Not described here, the proxy and environment construction (resp. env and proxy ) depend on the possible number of consecutive missed events (and rate of variation for event durations). However it is a simple event parse, creating an edge (in the trivial case) for each event of the mixed score.

top Groups

The IR resulted G is a set of parallelized machine G i describing the behavior of an actions sequence. So i is equal to the number of actions sequences of the mixed score. As env and proxy , evt parse the events of the mixed score. But it creates a machine G T i (called top group top ) for each events's no empty s.

TOP

where ev = evt(e, d, ) for the first rule and ev = evt(e, d, sa) for the second one. Is described here a step of the function ( evt ), note the step by passing ev from right to left of the parsed score (sc::sc'). A function ( top ), creating the related machine of sa (G sa ), is called (if s is not empty). The result is parallelized with the rest of the future groups G (built by the recursive call).

The second above rule assumes that we have associated a FSM G sa to a group sa triggered by the event e, using the predicate top . Implicitly, the group sa is treated like act(0, sa, [top]) where top is an adhoc strategy of sync and error handling specific to top-level groups.

Input 0 is the initialization of the accumulator, ev ::sc the rest of the score from the current event (included) and ẽ = closest(acc + d, sc) the next event after the event trigger. sa is the actions sequence parsed by the function top . The function closest associates to a duration d and a score sc the latest event e ∈ I in sc before d in the timeline defined by sc.

The trigger T τ associated to the signal or input event τ ∈ Σ sig ∪ I, and the finisher F i are defined as follows (for i = 2).

T τ = 0 1 2 τ ? τ ? and F 2 = 1 2
The finisher, set all seeker locations to stop locations. i can vary from 2 to 5 (@tight @global for a killable group).

al Groups

We will define here how a FSM is constructed during a parse of an actions sequence. The construction depends on the context, i.e the synchronization and error management attributes (al ) but a generic behavior can be defined. The current version of Antescofo defines the specific group top as a @loose @global group. A step is done by the construction of the FSM-part ( G), the accumulator incrementation (acc + d) and the update of the closest event (ẽ). A group parsing add a function ( al ), corresponding to the parsed group attributes, to construct its proper FSM (G sa ), and parallelizes it once terminated. Note that the function is initialized with the current accumulator updated, the rest of the mixed-score and the closest event.

We define the construction of the FSM for each attribute.

@loose @local

The loose local preserves the music coherence by following the musician only with the current tempo (delay in beats). It can be relevant for music accompaniment when it leads the part. Loose local groups have not late fire condition, if its trigger event is missed the group is not launch.

G(τ, d, [loose, local], sc, e) = 1 i 1 d τ !
To terminate a @loose @local group, F 2 or F 3 for killable groups, is used.

RR n°8701

2.8 @loose @global

The loose global behavior is quite as the loose local one, but for missed case the detected event (latter than the triggered one) current e split the group in two part. The past actions (from the trigger to the event) are directly fired, and other are played as a standard loose local group (= top group). The late loose global begin by a list of late actions and check if the next event is missed or not.

G(τ, d, [loose, global], sc, e) = 1 2 i 1 d τ ! 2 τ !
When e = closest(acc + d, sc), and otherwise, let e = closest(acc + d, sc),

G(τ, d, [loose, global], sc, e) = 1 2 i i ˜ i d e ? δ e ? 1 2 τ ! τ !
where δ is the difference between acc + d and the date of e in sc. To terminate a @loose @global group, F 2 or F 3 for killable groups, is used.

@tight @local

These attributes strictly synchronize the musician events played (during performance) with actions triggers. Tight machines manage standard and missed cases. We need one index more for the earlier events case related to the location l e . To terminate a @tight @local group, F 3 or F 4 for killable groups, is used.

@tight @global

The tight global machines regroup all the specifications of previous groups. We have a great view of what is the Antescofo's behavior on a specific case: from top to bottom a branch specify: standard, earlier standard, earlier missed and missed case. To terminate a @tight @global group, F 4 or F 5 for killable groups, is used.

special cases

Since actions can have a delay of 0, d and δ transitions can be erased (all parallel awaits (as e ?) are deleted too with their branches). An example is done for a @tight @global group.

RR n°8701

A jump can be seen for tight groups. It appends when an action delay is greater than a score event duration. 3 events are important for a jump since the previous delay wait on e, e * is waited to know the start of the next action's delay that can have the earlier event e played. where e * , e = closest(acc + d, sc).

For example for the score: we have e = e2, e * = e4, e = e5, when a1 is parsed.

Kill actions are managed as an action (the difference is the symbol that is an internal symbol (not an output)). The variation is on targeted groups (called killable) that have in parallel an asap transition waiting for the kill symbol. When fired the machine goes directly on a final state "kill".

Figure 1 :

 1 Figure 1: Grammar of the handled AST-mixed score

  all ∅ : A ∅ empty env sc : E proxy sc : P evt [ ], sc : G all sc : M all With M = E S and S = P || G

  evt sc::ev , sc : G evt sc, ev ::sc : G TOP evt sc::ev , sc : G 0, ev ::sc , ẽ top sa : G sa evt sc, ev ::sc : G T e + G sa + F i

  acc + d, sc, closest(acc + d, sc) al sa : G sa acc, sc, e al act(d, a, al )::sa : G(a, d, al , sc, e) + G sa msg acc + d, sc, ẽ al sa : G sa acc + d, sc, ẽ al sa : G sa acc, sc, e al act(d, sa , al )::sa : G + G sa T (g ) + G sa grp where |sa | > 1, g is a fresh signal name associated with the group sa , ẽ = closest(acc + d, sc), and G = G(g , d, al , sc, e) is a generic function that for a pair of attributes al constructs the corresponding part of FSM (detailed latter by cases).

G

  (τ, d, [tight, local], sc, e) When e = closest(acc + d, sc), and otherwise, let e = closest(acc + d, sc).RR n°8701G(τ, d, [tight, local], sc, e)

G

  (τ, d, [tight, global], sc, e) When e = closest(acc + d, sc), and otherwise, let e = closest(acc + d, sc).G(τ, d, [tight, global], sc, e)

G

  (τ, d, [tight, global], sc, e)

  BPM 144 evt(e 1 , 1, s 1 ); evt(e 2 , 0.5, []); evt(e 3 , 0.5, []); evt(e 4 , 2, []); evt(e 5 , 1, []); wheres 1 = act(0, s 2 , [tight; global]) s 2 = act(0.2, [a 0 ], [ ]); act(3, [a 1 ], [ ]);

RR n°8701