Charles-Gérard Lucas

Barbara Pascal
email: barbara.pascal@cnrs.fr

Nelly Pustelnik

Patrice Abry

C.-G Lucas

Hyperparameter selection for Discrete Mumford-Shah *

This work focuses on a parameter-free joint piecewise smooth image denoising and contour detection. Formulated as the minimization of a discrete Mumford-Shah functional and estimated via a theoretically grounded alternating minimization scheme, the bottleneck of such a variational approach lies in the need to fine-tune their hyperparameters, while not having access to ground truth data. To that aim, a Stein-like strategy providing optimal hyperparameters is designed, based on the minimization of an unbiased estimate of the quadratic risk. Efficient and automated minimization of the estimate of the risk crucially relies on an unbiased estimate of the gradient of the risk with respect to hyperparameters. Its practical implementation is performed using a forward differentiation of the alternating scheme minimizing the Mumford-Shah functional, requiring exact differentiation of the proximity operators involved. Intensive numerical experiments are performed on synthetic images with different geometry and noise levels, assessing the accuracy and the robustness of the proposed procedure. The resulting parameter-free piecewise-smooth estimation and contour detection procedure, not requiring prior image processing expertise nor annotated data, can then be applied to real-world images.

Introduction

Context -Image processing is characterized by several key tasks such as image recovery (e.g., debluring and/or denoising), feature extraction, segmentation, and contour detection, to name a few. To provide the user with the requested information, it is standard to perform successively a certain number of these tasks. A first major drawback of cascading tasks, is that important information might be thrown away at each stage. A second key issue is that each task might introduce estimation variance and/or regularization bias, which may accumulate and lead to subsequent errors on the target quantity. Finally, the selection of hyperparameters, e.g., regularization parameters, needs to be performed for each task independently, which might turn sub-optimal overall in minimizing the final error on the output estimate.

The benefit of performing jointly several steps has been illustrated in the context of texture segmentation [START_REF] Pascal | Strongly convex optimization for joint fractal feature estimation and texture segmentation[END_REF], providing a comparison between a two-step procedure (extract relevant local texture features followed by segmentation) against an original single-step procedure intertwining the estimation of relevant features and the segmentation procedure. Both strategies lead to strongly convex optimization schemes and fair comparisons can be provided by having recourse to an automatic hyperparameters selection procedure relying on Stein Unbiased Risk Estimator [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF].

Along this line, recent contributions in the image processing literature were dedicated to joint image denoising/restoration and contour detection [START_REF] Kiefer | PALMS Image Partitioning-A New Parallel Algorithm for the Piecewise Affine-Linear Mumford-Shah Model[END_REF][START_REF] Storath | Jump-sparse and sparse recovery using Potts functionals[END_REF][START_REF] Zach | Discretized convex relaxations for the piecewise smooth Mumford-Shah model[END_REF] via biconvex proximal optimization schemes. However, the automatic tuning of hyperparameters in this context has not been dealt with yet: this is the object of the present contribution. D-MS -This work focuses on a bi-convex formulation, refered as Discrete Mumford-Shah (D-MS) functional, that can trace back to the Mumford-Shah [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] or Geman and Geman functionals [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF], aiming to perform joint image denoising and contour detection, which may be written in the discrete variational formulation setting as:

min u∈R |Ω| ,e∈R |E| 1 2 u -z 2 2 + β (1 -e) Du 2 2 + λh(e), (1)
where z = u + σζ ∈ R |Ω| with ζ ∼ N (0 |Ω| , I |Ω| ) denotes the observed degraded image, defined on a grid of pixels Ω such that |Ω| = N , and σ > 0 is the known standard-deviation of the noise. The variable e is a discrete field defined on the lattice E, encapsulating the contour information, whose values are 1 when a contour is detected and 0 otherwise, D : R |Ω| → R |E| is a discrete difference operator such that Du lives on a lattice of contours E, denotes the component-wise product, h denotes a convex separable function enforcing sparsity having its minimum in 0 and such that, for every e = (e i ) 1≤i≤E , h(e) = i h i (e i ). β > 0 and λ > 0 are regularization parameters. The minimization is performed with SL-PAM, a nonconvex alternated minimization scheme, with descent parameters genuinely chosen to ensure fast convergence [START_REF] Foare | Semi-linearized proximal alternating minimization for a discrete Mumford-Shah model[END_REF] and recalled in Algorithm 4 in the Appendix 4. An exhaustive state-of-the-art regarding variations of the minimization problem [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF] is provided in Appendix D.

Hyperparameter selection -The aforementioned procedure for image denoising and contour detection involves hyperparameters, e.g. β and λ in [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF]. To reach satisfactory performance, the fine-tuning of these parameters is crucial. Although central in signal and image processing, this difficult task is still an ongoing challenge, particularly for variational methods.

The difficult problem of the selection of the regularization parameters of the D-MS functional is addressed considering a Stein Unbiased Risk Estimate (SURE), combined with a Finite Difference Monte Carlo (FDMC) strategy making its practical computation tractable. For a detailed state-of-the-art dedicated to hyperparameter selection, the reader could refer to Appendix C. The optimal regularization parameters obtained by minimizing FDMC SURE via exhaustive grid search are shown to lead to denoised estimates with high signal-to-noise ratio and relevant contours.

Then, to provide a fast procedure selecting the regularization parameters, a Stein Unbiased GrAdient Risk estimate (SUGAR) adapted to D-MS functional [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF] is designed, involving the Jacobian of the parametric estimator obtained from [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF]. Practical implementation of SUGAR requires iterative differentiation of the SL-PAM minimization scheme, for which closed-form formulas are provided. An averaging Monte Carlo strategy is discussed, providing a robust FDMC SUGAR estimator. The resulting procedure compares favorably against exhaustive grid search in terms of signal-to-noise ratio, while requiring a significantly smaller computational cost. To the best of our knowledge, the proposed automated D-MS bi-level scheme constitutes a first automated, prior-free and fast discrete Mumford-Shah-like formalism with automated selection of regularization parameters. Outline -The proposed automated and fast procedure is described in Section 2. Numerical experiments are provided in Section 3.

Hyperparameter selection for D-MS

Stein estimators for D-MS

The minimization of the D-MS functional of Eq. ( 19) provides both a piecewise smooth image reconstruction, denoted u(z; β, λ) and a set of detected contours, encapsulated into e(z; β, λ), depending on the choice of hyperparameters Θ = (β, λ) ∈ R + × R + . In such a context, we should ideally minimize a global error criterion

Θ ∈ Argmin Θ d( x(z; Θ), x), (2)
measuring the ability of the piecewise-smooth image and contour estimates x(z; Θ) = ( u(z; β, λ), e(z; β, λ)) to approximate the original data x = (u, e), using a measure of similarity d. If d is chosen to be a quadratic risk, it reads:

d( x(z; Θ), x) = E[ u(z; Θ) -u 2 2 ] + ζE[ e(z; Θ) -e 2 2 ], (3)
where ζ ≥ 0. However, the degradation model only describes how the observed image z relates to the ground truth image u, no prior knowledge about how the ground truth contours e are affected by the observation noise being assumed. Further, measuring the accuracy of the contours is a tedious task, involving complicated criteria, such as the Jaccard index [START_REF] Jaccard | Distribution de la flore alpine dans le bassin des Dranses et dans quelques régions voisines[END_REF]. For these reasons, in the present work, the quadratic error on which the choice of hyperparameter relies is chosen to be the quadratic estimation error on the reconstructed image (i.e. ζ = 0). The present work focuses on a strategy combining Finite Difference approximated differentiation and Monte Carlo averaging, which was first described by [START_REF] Ramani | Monte-carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms[END_REF]. Making use of a Finite Difference step > 0 and a Monte Carlo vector δ ∈ R N drawn from N (0 N , I N ), Finite Difference Monte Carlo (FDMC) SURE is defined as:

SURE ,δ (z; Θ|σ 2 ) := ( u(z; Θ) -z) 2 2 + 2 u(z + δ; Θ) -u(z; Θ), σ 2 δ -σ 2 N. (4)
It involves the denoised image u(z; β, λ), obtained from the minimization of [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF], and the standard deviation of the noise σ. Under technical assumptions detailed in Appendix E, the true inaccessible quadratic risk estimator (3) when ζ = 0 satisfies the following asymptotic unbiasedness property: lim

-→0 E[SURE ,δ (z; Θ|σ 2 )] = E[ u(z; Θ) -u 2 2 ]. (5)
Then, the design of a fast gradient-based hyperparameter selection strategy providing optimal hyperparameter from the minimization of ( 25) requires an unbiased estimate ∂ Θ SURE ,δ (z; Θ|σ 2 ) of the gradient of the quadratic risk with respect to hyperparameters Θ. Such a general procedure is sketched in Algorithm 5 and relies on the FDMC Stein Unbiased GrAdient Risk (SUGAR) estimate defined as:

SUGAR ,δ (z; Θ|σ 2 ) = 2∂ Θ u(z; Θ) * ( u(z; Θ) -z) + 2 (∂ Θ u(z + δ; Θ) -∂ Θ u(z; Θ)) * σ 2 δ, (6)
where ∂ Θ u(z; Θ) denotes the Jacobian of the parametric estimator u(z; Θ) with respect to the hyperparameters Θ. The first proposal of such Stein Unbiased GrAdient Risk (SUGAR) estimate was formulated by [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF] for i.i.d. Gaussian noise, and then extended in [START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF] for correlated noise. The main difficulty when it comes to practical implementation is to evaluate the Jacobian matrices.

Algorithm 1 Automated selection of hyperparameters. Input: Data z, , δ, and true or estimated σ 2 . Initialization:

Set Θ [0] ∈ R L . For t = 0 to T max -1 do       Compute SURE ,δ (z; Θ [t] |σ 2 ) Compute SUGAR ,δ (z; Θ [t] |σ 2 )
Update Θ [t] to Θ [t+1] via a gradient descent step

Output: Θ * = Θ [Tmax]

Differentiated SL-PAM

Practical evaluation of the risk and gradient of the risk estimates from Eq. ( 25) and ( 6), requires to compute the Jacobian of the D-MS estimator. No closed-form expression being available for u(z; β, λ), the derivatives are obtained from the iterative differentiation of the recursive scheme of DMS-SLPAM, Algorithm 4 in Appendix D. This strategy raises several technical issues. Indeed, following [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF] and [START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF], the Jacobian matrices ∂ Θ u(z; β, λ) and ∂ Θ u(z + δ; β, λ) are computed iteratively from a differentiated recursive scheme. Particularized to the case of D-MS estimates, the chain differentiation of the SL-PAM scheme (cf. Algorithm 4 in Appendix D) is derived in Algorithm 2. For ease of computation, a specific choice of the step-size d k = ηβ D 2 involved in the update of the variable e is considered, without inducing any loss of generality.

Algorithm 2 Iterative differentiation of SL-PAM

Input: Data z = {z, z + δ}. Set Θ = (β, λ) ∈ R + × R + . Initialization: u [0] = z, e [0] = 1 |E| , ∂ Θ u [0] = ∂ Θ u [0] = 0 N , ∂ Θ e [0] = ∂ Θ e [0] = 0 |E| . Set γ > 1 and η > 0. While |Ψ(u [k+1] , e [k+1] ) -Ψ(u [k] , e [k] )| > ξ                            Set c k = γβ D 2 and d k = ηβ D 2 u [k+1] = u (k) -1 c k ∇ u g(u [k] , e [k] ) u [k+1] = prox 1 c k f (•; z) u [k+1] Compute ∂ Θ u [k+1] from Eq. (8) Compute ∂ Θ u [k+1] from Eq. (9) For all i ∈ {1, . . . , |E|}              e [k] i = β(D i u [k+1] ) 2 + d k e [k] i 2 β(D i u [k+1] ) 2 + d k 2 e [k+1] i = prox λ 2β(D i u [k+1] ) 2 +d k hi ( e [k] i ) Compute ∂ Θ e [k+1] i from Eq. (10) Compute ∂ Θ e [k+1] i from Eq. (11)
General procedure -The purpose is to differentiate the mapping Θ → ( u(z; Θ), e(z; Θ)), where the estimates ( u(z; Θ), e(z; Θ)) are obtained solving [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF] for fixed z.

The recursive chain differentiation consists in differentiating step by step DMS-SLPAM each update of which can be written as v(z; Θ) = Γ(u(z; Θ), e(z; Θ), τ (Θ)), where u : R N ×R L → R N , e : R N × R L → R |E| and τ : R L → R are functions of the observed noisy image z and of the hyperparameters Θ, with respect to which the differentiation is to be performed and v(z; Θ) ∈ K, where K = R N when updating u or u, and K = R |E| when updating e or e.

Then, applying the chain rule differentiation principle yields the following partial derivative expression, for every component θ of the hyperparameter vector Θ and for every index j ∈ {1, . . . , dim(K)}:

∂ θ v j = N =1 (∂ u Γ j )(∂ θ u ) + |E| m=1 (∂ em Γ j )(∂ θ e m ) + (∂ τ Γ j )(∂ θ τ ) (7)
leading to the following closed form expression for

∂ Θ u [k+1] , ∂ Θ u [k+1] , ∂ Θ e [k+1] i and ∂ Θ e [k+1] i , for i ∈ {1, . . . , |E|}.
Iterative differentiation of DMS-SLPAM -Applying Formula (7) to each step of DMS-SLPAM leads to Algorithm 2. Proposition 1 provides closed-form expressions of the updates of the Jacobian matrices of the iterates involved in the minimization of a D-MS functional with h i = | • |, thus allowing an easy and direct implementation of FDMC SURE and FDMC SUGAR estimates of Eq. ( 25) and ( 6).

∂ θ u [k] = ∂ θ u [k] - 2β c k |E| i=1 (1 -e [k] i ) 2 D * i D i ∂ θ u [k] + 4β c k |E| i=1 (1 -e [k] i )∂ θ e [k] i D * i D i u [k] , (8)
∂ θ u [k+1] = c k c k + 1 ∂ θ u [k] + u [k] -z (βc + 1) 2 ∂ θ c k , (9)
where ∂ β c k = γ D 2 and ∂ λ c k = 0, and for every i ∈ {1, . . . , |E|}:

∂ θ e [k] i = 2D i u [k+1] D i ∂ θ u [k+1] d 2 (1 -e [k] i ) D i u [k+1] 2 + d 2 2 + d 2 ∂ θ e [k] i D i u [k+1] 2 + d 2 , (10)
∂ θ e [k+1] i = -∂ u φ [k+1] i ∂ θ u [k+1] e [k] i | e [k] i | I | e [k] i |>φ [k+1] i + ∂ θ e [k] i I | e [k] i |>φ [k+1] i - ∂ τ φ [k+1] i ∂ θ τ | e [k] i | e [k] i I | e [k] i |>φ [k+1] i , (11)
where

         ∂ u φ [k+1] i ∂ θ u [k+1] = -4τ Diu [k+1] Di∂ θ u [k+1] 2(Diu [k+1] ) 2 +d 2 , ∂ τ φ [k+1] i = 1 2(Diu [k+1] ) 2 +d , ∂ β τ = -λ β 2 , ∂ λ τ = 1 β . (12)
Proof. The proof is given in Appendix F and the notation I is defined in Appendix A.

Monte Carlo averaging strategy

Following [START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF], the risk and gradient of the risk FDMC Stein estimators, introduced in Eq. ( 25) and ( 6), are defined from one realization of the Monte Carlo vector δ. Yet, in the context of a parametric estimator u(z; Θ) obtained from the minimization of a nonconvex objective functional, such as [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF], it can be necessary to go further, and to consider Monte Carlo averaging strategies to get more robust risk and gradient of the risk estimates.

The Monte Carlo averaging strategy consists in averaging the FDMC Stein estimators of Eq. ( 25) and ( 6) over a certain number R of random Monte Carlo vectors δ (r) ∈ R N , independently sampled from the standard Gaussian distribution as stated properly in Definition 1.

Definition 1 (Monte Carlo averaged Stein estimators). Let z = u + σζ ∈ R |Ω| with ζ ∼ N (0 |Ω| , I |Ω| )
and let u(z; Θ) a parametric estimator of the underlying ground truth u, depending on some hyperparameters stored in Θ ∈ R L . For > 0 a Finite Difference step and ∆ = [δ (1) , . . . , δ (R) ] a concatenation of independent Monte Carlo vectors sampled from the standard Gaussian distribution. The Monte Carlo averaged SURE is defined as

SURE R ,∆ (z; Θ) := 1 R R r=1 SURE ,δ (r) (z; Θ), (13)
where SURE ,δ (r) (z; Θ) is the FDMC SURE [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Similarly, the Monte Carlo averaged SUGAR estimator writes

SUGAR R ,∆ (z; Θ) := 1 R R r=1 SUGAR ,δ (r) (z; Θ), (14)
involving SUGAR ,δ (r) (z; Θ), the FDMC SUGAR estimate (6).

-→0 E[SURE R ,∆ (z; Θ|σ 2 )] = Q[ u](Θ) (15)
and lim

-→0 E[SUGAR R ,∆ (z; Θ|σ 2 )] = ∂ Θ Q[ u](Θ). ( 16
)
Moreover, SUGAR R ,∆ (z; Θ|σ 2 ) is exactly the gradient of SURE R ,∆ (z; Θ|σ 2 ) with respect to the hyperparameters Θ.
Proof. The proof is provided in Appendix G.

Averaged SUGAR D-MS

The framework presented in Section 2, combined with Proposition 2, enable us to design an automated strategy to select the D-MS regularization parameters described below and assessed in Section 3. First, R independent Monte Carlo vectors δ (r) are sampled. The set ∆ = {δ (1) , . . . , δ (R) } is kept fixed throughout the procedure. Then, Algorithm 5 with the averaged estimates

Q(z; Θ|σ 2 ) = SURE R ,∆ (z; Θ|σ 2 ) and ∂ Θ Q(z, Θ|σ 2 ) = SUGAR R ,∆ (z; Θ|σ 2 )
, defined respectively in Eq. ( 13) and ( 14), whose practical implementation is based on Algorithm 2, provides the optimal hyperparameters. The overall procedure is referred to as Averaged SUGAR D-MS. Note that, for R = 1, ∆ = {δ (1) }, and one retrieves the standard SURE and SUGAR estimates. In the case when R = 1, the automated hyperparameter strategy is hence referred to as Standard SUGAR D-MS.

3 Performance assessment

Settings

Data -To assess the relevance of SURE [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] in the context of interface detection, as well as the efficiency of the proposed automated minimization making use of the SUGAR proposed in Section 2, systematic experiments are performed on the test data displayed in Fig. 7 and several noise levels are explored, corresponding to σ ∈ {0.01, 0.05, 0.1}. Additional experiments with other geometries and level of noise are provided in Appendix H. Algorithmic setup -See Supplementary materials. Performance criteria -In practice, standard and averaged SURE are compared to the following quadratic error :

Q( u|u) := u -u 2 2 .
To assess the performance of D-MS denoising with automatically selected hyperparameters, the quality of the reconstruction is quantified by the peak signal-to-noise ratio defined as: PSNR( u|u) = 20log 10 u u-u .

SURE for D-MS

We first illustrate in Fig. 2 the asymptotic unbiasedness of the standard and averaged SURE on the example z displayed in Fig. 7 (top-middle) with noise level σ = 0.05. To better locate and compare the minima, three level sets of SURE are displayed by the Matlab function contour.

Even though the overall shape of the standard SURE maps are similar to the quadratic error profile, Fig. 2(a-c) shows that the location of the minimum varies significantly with the Monte Carlo vector δ (r) . Averaged SURE also well reproduces the quadratic error map while being more robust to achieve the minimum (cf. Fig. 2(d)).

This first set of experiments illustrate that the proposed averaged SURE reproduces accurately the quadratic risk as expected from Proposition 2

Comparison between Standard and Averaged SUGAR D-MS

Fig. 3 investigates the ability of the hyperparameter selection strategies proposed in Section 2.4 for different numbers R ∈ {5, 10, 20} of Monte Carlo vector δ (r) to achieve the optimal hyperparameters minimizing the quadratic error Q( u|u).

The optimal hyperparameters Θ * (r) = (β * (r) , λ * (r) ) reached by the Standard SUGAR D-MS are scattered (Fig. 3 left), probably due to a lack of accuracy of the estimator SUGAR ,δ = SUGAR R=1 ,∆ . A first approach to alleviate the variability of the result is to carry out an averaging of R hyperparameters (β * (r) , λ * (r) ) obtained by the Standard SUGAR D-MS method: realizations are sufficient to achieve a good estimation of the optimal hyperparameters, second, we note that the proposed automated procedure is 20 times faster compared to exhaustive search, a grid search on averaged SURE requiring 60 minutes of calculation, while Averaged SUGAR D-MS requires 3 minutes, when using Matlab R2018a and an Intel Core i5 processor. Few estimated images obtained with the fully unsupervised parameter-free Averaged SUGAR D-MS are provided in Fig. 7 (2nd row).

Θ * R = (β * R , λ * R ) = 1 R R r=1 (β * (r) , λ * (r) ). (17)
SURE ,δ (1) (z; Θ|σ 2 ) SURE ,δ (2) (z; Θ|σ 2 ) SURE ,δ (3) (z; Θ|σ 2 ) SURE ,∆ (z; Θ|σ 2 ) E[ u(z; Θ) -u 2 2 ]
(d) (e) (f) (g) (h) (i) (j)

Real-world images

The proposed automated joint denoising and contour detection procedure Averaged SUGAR D-MS is evaluated on real-world images extracted from BSD69 dataset [START_REF] Roth | Fields of experts[END_REF] degraded with a Gaussian noise with σ = 0.05. In our experiments we set R = 5 and σ has been estimated from noisy data following Eq. ( 53) in Appendix H.3. Denoised images and contours provided by the proposed data-driven Averaged SUGAR D-MS strategy are compared with those yield by SUGAR T-ROF (a two-step procedure, consisting in, first, a piecewise constant denoising with automated tuning of the regularization parameter [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF], followed by an iterative thresholding procedure [START_REF] Cai | Multiclass segmentation by iterated ROF thresholding[END_REF]). In Fig. 8, we can observe that the denoising performance are very close for both procedures (in terms of SSIM, SUGAR T-ROF is slightly better) while the contour detection is significantly improved with the SUGAR D-MS procedure (which is confirmed when computing Jaccard index w.r.t contours obtained from the original image).

Conclusion

This work devises a procedure to automatically select the hyperparameters of the D-MS functional allowing to perform simultaneously image denoising and contour detection. This approach is fully unsupervised compared to alternative deep learning strategies such as [START_REF] Bertasius | Deepedge: A multi-scale bifurcated deep network for top-down contour detection[END_REF] and reference therein. However, in a future work, it would probably benefit to combine D-MS functional with unfolded deep learning strategies in order to design efficient supervised combined denoising and contour detection approaches.

A Matlab toolbox implementing the proposed automated image denoising and contour detection procedure is publicly available1 .

A Notations

Let H a real Hilbert space, and f : H → (+∞, +∞] a function which is proper, convex, and lower-semicontinuous and τ > 0 a real parameter, the proximity operator of τ f at point v ∈ H is uniquely defined by prox τ f (v) = arg min

u∈H 1 2 u -v 2 2 + τ f (v).
Additionally, let G be a real Hilbert space and let A : H → G a Lipschitzian map, we denote by L A > 0 the Lipschitz modulus of A, such that, for every (x, y) ∈ H×H, A(x)-A(y) ≤ L A x-y . Further, for every (x, y) ∈ R×R, we denote I x>y = 1 if x > y and 0 otherwise. Finally, I N denotes the identity matrix acting on R N , and 1 N (resp. 0 N ) is the vector of R N containing only ones (resp. zeros).

B State-of-the-art for contour detection in image processing

This work focuses on performing jointly piecewise smooth denoising and contour detection on images. In many classical approaches, image reconstruction is embedded into a variational formalism [START_REF] Golub | Tikhonov regularization and total least squares[END_REF][START_REF] Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF], which amounts to find a minimizer of a functional consisting of the sum of a data fidelity term and a prior penalization, i.e.,

minimize u 1 2 u -z 2 2 + γp(Du) (18)
where γ > 0, z ∈ R |Ω| denotes the observed degraded image, defined on a grid of pixels Ω, and D : R |Ω| → R |E| is a discrete difference operator such that Du lives on a lattice of contours E. Appropriate choices of the penalization term p, yield e.g. the Potts functional, when p = • 0 , or the Blake and Zisserman functional [START_REF] Blake | Visual reconstruction[END_REF][START_REF] Chambolle | Image segmentation by variational methods: Mumford and Shah functional and the discrete approximations[END_REF], corresponding to p(Du) = b min{ D b u q q , χ q }, for some q ∈ [1, ∞) and χ > 0, with D b being associated with several rows of D. In the same vein, considering a convex relaxation of Potts functional, contour detection can be obtained from the minimization of the Rudin-Osher-Fatemi (ROF) functional [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], which favors piecewise constant estimate when considering p(Du) = b D b u 2 . An alternative solution relies on a bi-convex formulation that can trace back to the Mumford-Shah [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] or Geman and Geman functionals [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF], which may be written in the discrete variational formulation setting as:

minimize u∈R |Ω| ,e∈R |E| 1 2 u -z 2 2 + β (1 -e) Du 2 2 + λh(e), (19)
where denotes the component-wise product, h denotes a convex function enforcing sparsity and β > 0 and λ > 0 are regularization parameters. This Discrete Mumford-Shah (D-MS) functional provides a piecewise-smooth reconstructed image u as well as estimated sparse contours e.

To achieve segmentation into K regions, Cai and Steidl designed an iterated thresolding strategy [START_REF] Cai | Multiclass segmentation by iterated ROF thresholding[END_REF] applied as a post-processing onto the minimizer of ROF functional. The resulting state-of-the-art two-step procedure, referred to as Thresholded ROF (T-ROF), was proven to be equivalent to minimizing the K-region piecewise constant Mumford-Shah functional. From this thresholded solution, it is then straightforward to identify the contours of the image. However, such an indirect contour extraction procedure restricts to closed contours. Fig. 5 shows a comparison between D-MS and T-ROF methods on a piecewise smooth image. The Mumford-Shah estimate is piecewise smooth preserving the discontinuities of the image while the ROF estimate is piecewise constant, leading to staircasing effects. We observe that T-ROF erroneously detects interfaces in areas on which the image is piecewise smooth, as opposed to the D-MS whose estimated contour variable is approximately zero everywhere except at the location of the actual signal discontinuity.

C State-of-the-art for Hyperparameter selection

All aforementioned procedures for image denoising and contour detection involve hyperparameters, e.g. β and λ in [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF]. To reach satisfactory performance, the fine-tuning of these parameters is crucial. Although central in signal and image processing, this difficult task is still an ongoing challenge, particularly for variational methods.

A first class of methods relying on hierarchical Bayesian approaches and has been widely used, both in signal and image processing [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF][START_REF] Dobigeon | Joint segmentation of piecewise constant autoregressive processes by using a hierarchical model and a Bayesian sampling approach[END_REF][START_REF] Molina | Image restoration in astronomy: A Bayesian perspective[END_REF][START_REF] Vacar | Unsupervised joint deconvolution and segmentation method for textured images: a Bayesian approach and an advanced sampling algorithm[END_REF]. The drawbacks of Bayesian methods are that they rapidly become computationally heavy as the model for observed data gets more complicated, and their computational cost increases with the number of hyperparameters to be tuned. For specific 1D denoising problems, efficient hybrid variational/Bayesian strategies can be designed [START_REF] Frecon | Bayesian Selection for the 2 -Potts Model Regularization Parameter: Constant Signal Denoising[END_REF].

Several other classes of methods, such as cross-validation or Stein Unbiased Risk Estimate (SURE) formulation, can be formulated as a bilevel optimization problem. Cross-validation relies on a given labeled data set composed of noisy samples with their associated ground truth [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF][START_REF] Stone | Cross-validation: A review[END_REF]. However, in several real-world applications, such as medical imaging [START_REF] Marin | Mammographic evidence of microenvironment changes in tumorous breasts[END_REF] or nonlinear physics problems [START_REF] Pascal | Parameter-free and fast nonlinear piecewise filtering: application to experimental physics[END_REF], obtaining a large enough labeled dataset is very challenging, if not impossible. Hence, SURE, initially proposed in [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF], has long been favored for its combined simplicity and efficiency. Stein-based hyperparameter strategies rely on an additive Gaussian noise model to design an estimate of the inaccessible true risk, defined as the quadratic error between the estimate and ground truth. The major advantage of these approaches is that they do not require to access ground truth. Then, the selection of optimal hyperparameters is done by minimizing SURE and by making use of Finite Difference strategies [START_REF] Shen | Adaptive model selection[END_REF][START_REF] Ye | On measuring and correcting the effects of data mining and model selection[END_REF] or/and Monte Carlo averaging [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF][START_REF] Ramani | Monte-carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms[END_REF][START_REF] Girard | A fast 'Monte-Carlo cross-validation' procedure for large least squares problems with noisy data[END_REF], to yield tractable and fast implementation of Stein-based risk estimates.

However, the strategy to find the optimal hyperparameters for a specific criterion has a huge impact on the solution both in terms of quality assessment and in terms of computational load. The most standard approach consists in computing a chosen error criterion over a grid of parameters [START_REF] Ramani | Monte-carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms[END_REF][START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF], and to select the parameter of the grid for which the error is minimal. Such a grid search procedure suffers from a high computation cost, especially when dealing with L ≥ 2 regularization parameters. To circumvent this difficulty, efficient automated minimization methods are required. It was early envisionned by Chaux et al. [START_REF] Chaux | A nonlinear Stein-based estimator for multichannel image denoising[END_REF], who proposed and assessed numerically an empirical descent algorithm for automated choice of regularization parameters, but with no convergence guarantee. A deeper theoretical analysis was then provided by Deledalle et al. [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF], evidencing sufficient conditions so that Stein Unbiased Risk Estimate is differentiable with respect to hyperparameters, thus enabling to define the Stein Unbiased GrAdient of the Risk (SUGAR) estimator and to provide a practical implementation based on an iterative differentiation strategy. Combining SUGAR with a quasi-Newton descent procedure, a fast algorithm was designed to achieve optimal hyperparameters selection for objective functions of the form [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF]. This strategy, later extended in [START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF][START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF] for correlated noise, proved its efficiency for texture segmentation [START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF], piecewise linear signal denoising [START_REF] Pascal | Parameter-free and fast nonlinear piecewise filtering: application to experimental physics[END_REF], and in spatial-spectral deconvolution for large multispectral data [START_REF] Ammanouil | A parallel and automatically tuned algorithm for multispectral image deconvolution[END_REF].

D Minimization of the discrete Mumford-Shah functional

The D-MS functional introduced in Eq. ( 19) being nonconvex, standard proximal algorithms [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Parikh | Proximal algorithms[END_REF] cannot be used directly for its minimization. However, the fact that the functional is separately convex with respect to each variable advocates the use of alternating schemes. Among the vast variety of existing alternating algorithms benefiting from convergence guarantees [START_REF] Foare | Semi-linearized proximal alternating minimization for a discrete Mumford-Shah model[END_REF][START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], a numerically efficient procedure for the minimization of D-MS like functionals appears to be the Semi-Linearized Proximal Alternating Minimization (SL-PAM) scheme proposed in [START_REF] Foare | Semi-linearized proximal alternating minimization for a discrete Mumford-Shah model[END_REF], whose iterations in the general setting of Problem 1 are recalled in Algorithm 3. (

) 20
The algorithmic scheme SL-PAM (Algorithm 3) is an hybrid version between PAM [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF] and PALM [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. The key ingredient for the efficiency of SL-PAM consists in avoiding the linearization with respect to the variable e [k] , enabling to choose larger descent steps. Under some technical assumptions, such as the existence of a closed-form expressions of the involved proximity operators, the sequence u [k] , e [k] k∈N converges toward a critical point of Ψ(u, e).

Algorithm 3 SL-PAM

Initialization: u [0] = z, e [0] = 1 |E| , γ > 1 and ξ > 0. While |Ψ(u [k+1] , e [k+1] ) -Ψ(u [k] , e [k] )| > ξ         Set c k = γL β∇ug(•,e [k] ) and d k > 0 u [k] = u [k] -β c k ∇ u g(u [k] , e [k] ) u [k+1] = prox 1 c k f u [k] e [k+1] = prox 1 d k (λh+βg(u [k+1] ,•)) e [k]
The piecewise smooth image denoising and contour detection strategy defined by ( 19) and on which this paper focuses corresponds to a particularization of Problem 1. The three terms of the objective function Ψ of Eq. ( 20) are particularized to

                   f (u) = 1 2 u -z 2 2 , g(u, e) = |E| i=1 (1 -e i ) 2 (D i u) 2 , h(e) = |E| i=1 h i (e i ) (21)
where, for all i ∈ {1, . . . , |E|}, D i denotes the i th -row of the discrete gradient operator D, and h i : R → (-∞, +∞] is a separable proper, lower semi-continuous, and convex function having a proximal operator with known closed-form expression. The iterations of Algorithm 3 specified to the minimization of D-MS lead to Algorithm 4 [START_REF] Foare | Semi-linearized proximal alternating minimization for a discrete Mumford-Shah model[END_REF].

For a detailed discussion of the convergence behavior depending on the choice of the descent steps γ and d k , the reader is referred to [START_REF] Foare | Semi-linearized proximal alternating minimization for a discrete Mumford-Shah model[END_REF]. The most efficient setting appears to choose both of them the smallest possible.

E Risk estimation

As previously discussed in introduction, many variational approaches for image restoration and contour detection consists in designing a parametric estimator u(z; Θ), e.g., defined as a minimizer of ( 18) or [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF], which aims at providing the best possible estimate of a quantity of interest u from noisy observations z. By construction, the quality of this estimate crucially relies on the precise selection of the hyperparameters Θ, which can be for instance the regularization parameters β and λ in D-MS functional [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF].

Quadratic risk based parameter selection

The hyperparameters tuning task is commonly formulated as the minimization of the following quadratic risk :

Q[ u](Θ) = E[ u(z; Θ) -u 2 2 ], (22)
Algorithm 4 DMS-SLPAM to solve [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF] Input: Data z. Set β > 0, λ > 0. Initialization:

u [0] = z, e [0] = 1 |E| ∈ R |E| . Set γ > 1 and ξ > 0. While |Ψ(u [k+1] , e [k+1] ) -Ψ(u [k] , e [k] )| > ξ                   Set c k = γβ D 2 and d k > 0. u [k] = u [k] -β c k ∇ u g(u [k] , e [k] ) u [k+1] = prox 1 c k f ( u [k] ) For all i ∈ {1, . . . , |E|}          e [k] i = β(D i u [k+1] ) 2 + d k e [k] i 2 β(D i u [k+1] ) 2 + d k 2 e [k+1] i = prox λ 2β(D i u [k+1] ) 2 +d k hi ( e [k] i )
measuring the expected reconstruction error made when estimating ground truth u by u(z; Θ). The expectation in Eq. ( 22) runs over the realizations of the noise corrupting z.

In practice, u being unknown and the number of observed samples z being limited, if not reduced to one, the exact quadratic risk Q[ u](Θ) of Eq. ( 22) is not accessible. Thus, the minimization of the quadratic risk Q[ u](Θ) is replaced by the minimization of some estimate Q(z; Θ|σ 2 ) computed from a single noisy sample z, not requiring the knowledge of ground truth but only some prior knowledge about the noise, e.g., its standard deviation σ :

Θ ∈ Argmin Θ Q(z; Θ|σ 2 ). (23)
Then, the design of a fast gradient-based hyperparameter selection strategy providing optimal hyperparameters from the minimization of ( 23) requires an unbiased estimate ∂ Θ Q(z; Θ|σ 2 ) of the gradient of the quadratic risk with respect to hyperparameters Θ. Such a general procedure is sketched in Algorithm 5.

Algorithm 5 Automated selection of hyperparameters. Input: Data z and true or estimated σ 2 . Initialization:

Set Θ [0] ∈ R L . For t = 0 to T max -1 do       Compute Q(z; Θ [t] |σ 2 ) Compute ∂ Θ Q(z; Θ [t] |σ 2 ) Update Θ [t] to Θ [t+1] via a gradient descent step Output: Θ * = Θ [Tmax]
Stein Unbiased Risk Estimate -To address the fact that the ground truth u is unknown, the pioneer work of Stein [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] proposed an unbiased estimate of the quadratic risk, based on an i.i.d. Gaussian noise additive model in which the observations are supposed to write

z = u + σζ, ζ ∼ N (0 N , I N ) (24)
with N = |Ω| is the number of pixels and σ 2 the known variance of the noise. Then, under integrability and regularity assumptions, together with the observation model ( 24), the so-called Stein Unbiased Risk Estimator (SURE) was derived in [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF], and has then been intensively used in signal and image processing [START_REF] Pascal | Strongly convex optimization for joint fractal feature estimation and texture segmentation[END_REF][START_REF] Ammanouil | ADA-PT: An adaptive parameter tuning strategy based on the weighted stein unbiased risk estimator[END_REF][START_REF] Chaux | A nonlinear Stein-based estimator for multichannel image denoising[END_REF][START_REF] Pascal | Parameter-free and fast nonlinear piecewise filtering: application to experimental physics[END_REF][START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF]. In most applications, the original Stein estimator is not usable directly and further strategies are necessary to yield a practical estimator. The present work focuses on a strategy combining Finite Difference approximated differentiation and Monte Carlo averaging, which was first described by [START_REF] Ramani | Monte-carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms[END_REF]. Making use of a Finite Difference step > 0 and a Monte Carlo vector δ ∈ R N drawn from N (0 N , I N ), Finite Difference Monte Carlo (FDMC) SURE is defined as:

SURE ,δ (z; Θ|σ 2 ) := ( u(z; Θ) -z) 2 2 + 2 u(z + δ; Θ) -u(z; Θ), σ 2 δ -σ 2 N, (25)
Under the Lipschitzianity with respect to z of u(z; Θ) and the natural unambiguity property u(0 N , Θ) = 0 N , the true inaccessible quadratic risk estimator [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] satisfies the following asymptotic unbiasedness property:

lim -→0 E[SURE ,δ (z; Θ|σ 2 )] = Q[ u](Θ), (26)
where the expectation is to be understood on both the realizations of the observation noise ζ appearing in Eq. ( 24), and the realizations of the Monte Carlo vector δ. Eq. ( 26) ensures that, for small enough Finite Difference step , and provided that N is large enough so that the Monte Carlo strategy is relevant, a minimizer of SURE ,δ (z; Θ|σ 2 ) is an approximately optimal set of hyperparameters in terms of quadratic risk.

Risk estimate minimization -The gradient-based strategy sketched at Algorithm 5 when

Q(z; Θ|σ 2 ) = SURE ,δ (z; Θ|σ 2 ) (27)
relies on the FDMC Stein Unbiased GrAdient Risk (SUGAR) estimate defined as:

SUGAR ,δ (z; Θ|σ 2 ) = 2∂ Θ u(z; Θ) * ( u(z; Θ) -z) + 2 (∂ Θ u(z + δ; Θ) -∂ Θ u(z; Θ)) * σ 2 δ, (28)
where ∂ Θ u(z; Θ) denotes the Jacobian of the parametric estimator u(z; Θ) with respect to the hyperparameters Θ. The first proposal of such Stein Unbiased GrAdient Risk (SUGAR) estimate was formulated by [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF] for i.i.d. Gaussian noise, and then extended in [START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF] for correlated noise. The main difficulty when it comes to practical implementation is to evaluate the Jacobian matrices.

In [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF][START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF], the authors proposed an efficient implementation when u(z; Θ) is estimated from the resolution of a convex minimization problem of the form [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF] while in this contribution we extend it in the context of interface detection involving a minimization problem such as [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF] solved with SL-PAM described in Section D.

Under technical assumptions such as Lipschitzianity of u(z; Θ) with respect to Θ and z, it has been proved in [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF] that the quadratic risk estimator ( 25) is weakly differentiable with respect to Θ and its gradient is exactly the gradient estimator recalled in [START_REF] Donoho | Denoising by soft-thresholding[END_REF], i.e.,

∂ Θ SURE ,δ (z; Θ|σ 2 ) = SUGAR ,δ (z; Θ|σ 2 ). (29)
Eq. ( 29) ensures that the gradient estimate SUGAR ,δ (z; Θ|σ 2 ) is indeed the gradient of the quadratic risk estimate SURE ,δ (z; Θ|σ 2 ) with respect to hyperparameters Θ, justifying the use of the gradient descent approach of Algorithm 5 to solve a particular instance of Problem [START_REF] Chambolle | Image segmentation by variational methods: Mumford and Shah functional and the discrete approximations[END_REF] when Q(z; Θ|σ 2 ) is defined by [START_REF] Dobigeon | Joint segmentation of piecewise constant autoregressive processes by using a hierarchical model and a Bayesian sampling approach[END_REF]. Additionally, FDMC SUGAR estimator introduced in ( 28) is an asymptotically unbiased estimator of the gradient of the true quadratic risk, i.e.

lim

-→0

E[SUGAR

,δ (z; Θ|σ 2 )] = ∂ Θ Q[ u](Θ), (30)
where

∂ Θ Q[ u](Θ)
is the true inaccessible gradient of quadratic risk with respect to hyperparameters Θ, and the expectation is to be understood on both the realizations of the observation noise ζ appearing in Eq. ( 24) and the realizations of the Monte Carlo vector δ. The asymptotic unbiasedness of the gradient estimate ensures that the risk profile around its minimum is well enough reproduced by Stein-like estimates so that Algorithm 5 can be reasonably supposed to output a good approximation of the true optimal hyperparameters.

F Iterative differentiation of SL-PAM for D-MS

F.1 Update of ∂ θ ũ[k]
For the function g given in Eq. ( 21), the update rule of u [k] reads:

u [k] = u [k] - β c k ∇ u g(u [k] , e [k] ) = u [k] - 2β c k |E| i=1 (1 -e [k] i ) 2 D * i D i u [k] . (31)
The update of ũ[k] can be written:

ũ[k] = Γ(u [k] , e [k] , τ [k] ), (32)
where

       Γ(u, e, τ ) = u -τ |E| i=1 (1 -e i ) 2 D * i D i u, τ [k] = 2β c k . (33)
The derivative ∂ θ v of v = Γ(u, e, τ ), for θ ∈ {β, λ}, is:

∂ θ v =∂ θ u -τ |E| i=1 (1 -e i ) 2 D * i D i ∂ θ u + 2τ |E| i=1 (1 -e i )∂ θ e i D * i D i u -∂ θ τ |E| i=1 (1 -e i ) 2 D * i D i u, (34)
and, since c k = βγ D 2 , for θ ∈ {β, λ},

∂ θ τ [k] = ∂ θ 2β c k = ∂ θ 2 γ D 2 = 0. (35)
F.2 Update of ∂ θ u [k+1] The function f given in Eq. ( 21) has a proximal operator with a closed form expression. Thus the update rule of u [k] can be explicitely expressed as follows:

u [k+1] = prox 1 c k f ( u [k] ) = c k u [k] + z c k + 1 . ( 36
)
For the update of u [k+1] , we thus have:

u [k+1] = Γ(ũ [k] , 0 |E| , τ [k] ), (37)
where

Γ(u, e, τ ) = τ u+z τ +1 , τ [k] = c k . ( 38
)
The derivative ∂ θ v of v = Γ(v, e, τ ), for θ ∈ {β, λ}, is:

∂ θ v = τ τ + 1 ∂ θ u + u -z (τ + 1) 2 ∂ θ τ, (39)
and

∂ β τ [k] = γ D 2 and ∂ λ τ [k] = 0. F.3 Update of ∂ θ ẽ[k]
In Algorithm 4, the parameter for the update of e [k+1] is set to d k = βd where d = η D 2 2 . This choice is discussed in [START_REF] Foare | Semi-linearized proximal alternating minimization for a discrete Mumford-Shah model[END_REF] and gives good numerical results for some values of η. This setting simplifies the computation of the derivatives due to the linear dependance of d k with β. Indeed, the update rule of e i , for every i ∈ {1, . . . , |E|}, can be rewritten as follows:

e [k] i = (D i u [k+1] ) 2 + de [k] i 2 (D i u [k+1] ) 2 + d 2 . (40)
Thus, for every i ∈ {1, . . . , |E|},

ẽ[k] i = Γ i (u [k+1] , e [k] , τ [k] ), (41)
where

Γ i (u, e, τ ) = (Diu) 2 + τ 2 ei (Diu) 2 + τ 2 , τ [k] = d. (42)
The derivative ∂ θ v of v = Γ(u, e, τ ), for θ ∈ {β, λ} and for i ∈ {1, . . . , |E|}, is:

∂ θ v i = 2D i uD i ∂ θ u (D i u) 2 + τ 2 - 2D i uD i ∂ θ u (D i u) 2 + τ 2 e i (D i u) 2 + τ 2 2 + τ 2 ∂ θ e i (D i u) 2 + τ 2 + ∂ θ τ 2 e i (D i u) 2 + τ 2 - [(D i u) 2 + τ 2 e i ] ∂ θ τ 2 (D i u) 2 + τ 2 2 , (43)
and ∂ θ τ For the update of e

[k+1] i

, for every i ∈ {1, . . . , |E|}, the fonction h in Eq. ( 21) is chosen with h i = |.|. This latter function corresponds to the common 1 -norm penalization of the contour. The setting d k = βd simplifies this update which now reads:

e [k+1] i = prox φ [k+1] i |.| ( e [k] i ), φ [k+1] i = φ i (u [k+1] ; τ ), (44) where 
  φ i (u, τ ) = τ [2(Diu) 2 +d] , τ = λ β , prox φi(u,τ )|.| (e i ) = max(0, 1 -φi(u;τ ) |ei| )e i . (45)
For every i ∈ {1, . . . , |E|}, e

[k+1] i = Γ i (u [k+1] , ẽ[k] , τ [k] ), (46)
where

Γ i (u, e, τ ) = prox φi(u;τ )|.| (e i ), τ [k] = λ β . (47)
The derivative ∂ θ v of v = Γ(u, e, τ ), for θ ∈ {β, λ} and for i ∈ {1, . . . , |E|}, is:

∂ θ v i = -∂ u φ i ∂ θ u e i |e i | I |ei|>φi(u,τ ) + ∂ θ e i I |ei|>φi(u,τ ) - ∂ θ τ 2 (D i u) 2 + d e i |e i | I |ei|>φi(u,τ ) , (48)
with Jacobian matrices product

∂ u φ i ∂ θ u = - τ 2 (D i u) 2 + d 2 (4D i uD i ∂ θ u) , (49)
and

∂ β τ [k] = 1 β and ∂ λ τ [k] = -λ β 2 .

G Proof of Proposition 2

For each Monte Carlo vector δ (r) , the FDMC SURE ,δ (r) and SUGAR ,δ (r) estimates, defined at Eq. ( 25) and ( 28

H Additional experiments

H.1 Algorithmic setup SL-PAM -The minimization of the D-MS functional [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF] providing estimates of both the piecewise smooth image and its salient contours, is performed running Algorithm 4. The stopping criterion, based on the objective function increments, is set to ξ = 10 -4 , while the descent steps Stein estimators -FDMC SURE ( 25) is computed with a Finite Difference step

= 2 σ N α , 0 < α < 1 ( 50
)
where σ is the standard deviation of the noise on the observed image z ∈ R N . Formula (50) derives from a heuristic reasoning developed in [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF], in the context of 1 -norm penalization. The dependency of the Finite Difference step on the size of the data is controlled via the exponent α, which is fixed at α = 0.3 for all the numerical simulations. In the systematic numerical experiments, four values of the number R of realizations of the Monte Carlo vector δ (r) are envisioned and systematically compared: R ∈ {1, 5, 10, 20}. BFGS algorithm -To perform the risk minimization described in Algorithm 5 for different choices of Q(z; Θ |σ2 ) and ∂ Θ Q(z; Θ |σ 2 ) we used the GRadient-based Algorithm for Non-Smooth Optimization, implemented in GRANSO toolbox 2 , consisting of the low memory BFGS quasi-Newton algorithm proposed in [START_REF] Curtis | A BFGS-SQP method for nonsmooth, nonconvex, constrained optimization and its evaluation using relative minimization profiles[END_REF] with box constraints, enabling to enforce positivity of β and λ. The maximal number of iterations of BFGS Algorithm 5 is set to T max = 20, while the stopping criterion on the gradient norm is set to 10 -8 . Further, it is well-documented that the initialization of quasi-Newton algorithms might drastically impact their convergence. Hence, we propose a model-based strategy for initializing Algorithm 5. Inspired from the initialization strategies proposed in [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF][START_REF] Pascal | Automated data-driven selection of the hyperparameters for Total-Variation based texture segmentation[END_REF], hyperparameters Θ = (β, λ) are initialized as

β (0) = N σ 2 4 Dz 2 2 and λ (0) = β (0) Dz 2 2 2N , (51)
while, for κ = 0.9, the initial approximated inverse Hessian involved in the BFGS strategy is set to

H (0) = diag κβ (0) ∂ β Q(z; Θ (0) |σ 2 ) , κλ (0)
∂ λ Q(z; Θ (0) |σ 2 ) . ( 52

H.2 Performance w.r.t noise level

We now focus on the Averaged SUGAR D-MS for R = 5 and assess its performance for the different geometries and noise levels displayed in Fig. 6. Averaged PSNR for 10 realizations of the noise are reported in Table 1, denoised images and detected contours are displayed in Fig. 7 (1st and 3rd rows) for one realization of the noise. As expected, the PSNR decreases as the noise level increases. Further, large error bars are observed mainly due to the variability of the gradient descent scheme in Algorithm 5, the realizations of the Monte Carlo vector, or the SL-PAM non-convex minimization procedure.

H.3 Impact of the estimation of σ

On real data, the noise level σ needed to implement Averaged SUGAR D-MS is unknown. The most usual method to estimate the noise standard deviation is the median absolute deviation (MAD) of 2D discrete wavelet coefficients [START_REF] Donoho | Denoising by soft-thresholding[END_REF]:

σ = MAD |ψ H,k |, |ψ V,k |, |ψ D,k ||k ∈ {1, . . . , N 4 } 0.6745 , (53)
where ψ H,k , ψ V,k , ψ D,k are the three wavelets coefficents (horizontal, vertical and diagonal) at the finest scale. Table 1 presents the PSNR values with the estimated noise level σ. In addition, some estimates are depicted in Fig. 7 (2nd and 4th rows).

The results with either estimated or true noise level are similar, attesting the robustness of Averaged SUGAR D-MS to the estimation of the noise level.

H.4 Real-world images

The proposed automated joint denoising and contour detection procedure Averaged SUGAR D-MS is evaluated on real-world images extracted from BSD69 dataset [START_REF] Roth | Fields of experts[END_REF] degraded with a Gaussian noise with σ = 0.05. In our experiments we set R = 5 and σ has been estimated from noisy data following (53). Denoised images and contours provided by the proposed datadriven Averaged SUGAR D-MS strategy are compared with those yield by SUGAR T-ROF (a two-step procedure, consisting in, first, a piecewise constant denoising with automated tuning of the regularization parameter [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF], followed by an iterative thresholding procedure [START_REF] Cai | Multiclass segmentation by iterated ROF thresholding[END_REF]). In Fig. 8, we can observe that the denoising performance are very close for both procedures (in terms of SSIM, SUGAR T-ROF is slightly better) while the contour detection is significantly improved with the SUGAR D-MS procedure (which is confirmed when computing Jaccard index w.r.t contours obtained from the original image).

Proposition 1 .

 1 Considering the D-MS functional (19) when h i = | • | and its minimization via SL-PAM Algorithm 3 in Appendix D with d k = βd, d = η D 2 2 , for every θ ∈ {β, λ}:

Proposition 2 .

 2 Let u(z; Θ) an estimator being uniformly Lipschitz w.r.t the observations z and w.r.t the hyperparameters Θ, with a Lipschitz modulus L u(z;•) independent of z, and satisfying the univocity condition u(0 N ; Θ) = 0 N . For a infinitesimal positive Finite Difference step and ∆ = δ (1) , . . . , δ (R) a collection of independent standard Gaussian Monte Carlo vectors, the Monte Carlo averaged estimates SURE R ,∆ (z; Θ) and SUGAR R ,∆ (z; Θ) are asymptotically unbiased estimates of respectively the risk and of the gradient of the risk with respect to hyperparameters lim

Figure 1 :

 1 Figure 1: Piecewise smooth grey level image corrupted by i.i.d. Gaussian noise with level σ ∈ {0.01, 0.05, 0.1}. Associated D-MS estimates u and e (superimposed in red) The D-MS hyperparameters are selected with the proposed Averaged SUGAR D-MS (with R = 5) using the true standard deviation σ.

Figure 2 :Figure 3 :

 23 Figure 2: Comparison between the quadratic error, standard and averaged SURE estimates for D-MS denoising of the image displayed in Fig. 7(middle). 1st row -Map on a logarithmic grid of 40 × 40 hyperparameters Θ = (β, λ): (a-c) SURE ,δ (r) (z; Θ|σ 2 ) values for some realizations of the Monte Carlo vector, (d) SURE ,∆ (z; Θ|σ 2 ) values for R = 5 realizations of the Monte Carlo vector and (e) quadratic error Q( u(z; Θ)|u) values with level sets (black lines). 2nd row -Optimal solutions u(z; Θ Grid ), e(z; Θ Grid ) obtained from a grid search over each map. The red (resp. yellow and green) cross corresponds to the solution displayed in (j) (resp. (f)-(h) and (i)) associated with the minimum of the quadratic error grid (e) (resp. SURE estimate grids (a)-(c) and (d)).

Figure 4 :

 4 Figure 4: Comparisons between SUGAR T-ROF and SUGAR D-MS for noisy images extracted from the BSD69 dataset [12].

Figure 5 :

 5 Figure 5: Comparison of state-of-the-art convex variational formulation T-ROF and the studied non-convex D-MS performing image denoising and contour extraction. From left to right: (a) Original noise-free piecewise smooth image, (b) Observations z corrupted by an additive Gaussian noise, (c) State-of-the-art ROF piecewise constant estimate and contours derived from thresholding into K = 3 regions (displayed in red), and (d) Studied D-MS piecewise smooth approximation and estimated contours (displayed in red).

Problem 1 (

 1 Nonconvex and nonsmooth minimization). Let f : R |Ω| → (-∞, +∞], h : R |E| → (-∞, +∞] two proper lower semi-continuous functions and g : R |Ω| × R |E| → (-∞ + ∞] a C 1 function. Let λ > 0 and β > 0. We aim to estimate: ( u, e) ∈ Argmin u∈R |Ω| ,e∈R |E| Ψ(u, e) := f (u) + βg(u, e) + λh(e).

F. 4

 4 [k] = 0, which yields to the result in[START_REF] Lucas | Hyperparameter selection for Discrete Mumford-Shah[END_REF] Proposition 2]. Update of ∂ θ e[k+1]

1 Figure 6 :

 16 Figure 6: Piecewise smooth grey level images (σ = 0) corrupted by i.i.d. Gaussian noise with level σ ∈ {0.01, 0.03, 0.05, 0.07, 0.1}.

  ) (a) σ = 0.01 (b) σ = 0.03 (c) σ = 0.05 (d) σ = 0.07 (e) σ = 0.1

Figure 7 :

 7 Figure 7: D-MS estimates u and e (superimposed in red) of piecewise smooth grey level images corrupted by i.i.d. Gaussian noise with noise level σ displayed in Fig. 6. The D-MS hyperparameters are selected with the proposed Averaged SUGAR D-MS using either the true standard deviation σ (first and third rows) or the estimated standard deviation σ (second and fourth rows).

  ) are asymptotically unbiased and SUGAR ,δ (r) , is the derivative of SURE ,δ (r) , w.r.t. Θ. Then, by linearity of both the limit lim →0 and the summation over the R Monte Carlo vectors, the Monte Carlo averaged estimates SURE

	R ,∆ (z; Θ|σ 2 ) and SUGAR	R ,∆ (z; Θ|σ 2 )
	are also unbiased and SUGAR R ,∆ (z; Θ|σ 2 ) is the derivative of SURE R ,∆ (z; Θ|σ 2 ) w.r.t. Θ.

Table 1 :

 1 PSNR values with 95% confidence interval for true noise level σ and for estimated noise level σ. .85 ± 0.06 43.64 ± 0.12 41.77 ± 0.06 41.08 ± 0.25 0.03 38.94 ± 0.12 38.89 ± 0.13 31.80 ± 2.73 33.64 ± 2.23 0.05 34.33 ± 0.70 34.71 ± 0.12 26.05 ± 2.92 26.37 ± 2.80 0.07 31.60 ± 0.52 30.55 ± 1.78 26.73 ± 2.46 22.64 ± 2.68 0.1 28.86 ± 0.47 27.24 ± 1.84 21.25 ± 3.03 22.61 ± 3.33

	σ	Losange True σ Estimated σ	True σ	Ellipse Estimated	σ
	0.01 43

https://github.com/charlesglucas/sugar_dms

http://www.timmitchell.com/software/GRANSO/

supported by ANR-19-CE48-0009s MULTISC-IN.


Christophe Maudoux
email: christophe.maudoux@cnam.fr

Selma Boumerdassi

LemonLDAP

Single Sign-On (SSO) is an authentication scheme that permits a user to log in with a single ID to any of several related, yet independent, systems. It allows the user to log in once and access services without re-entering authentication credentials. SSO solutions are classified depending on Authentication, Authorization, and Accounting features. AAA is an acronym defining a framework for intelligently controlling access to resources, enforcing security policies, auditing usage, and providing the information necessary to bill for services. These combined processes are considered important for effective network management and security.

Multi Factor Authentication (MFA) is also supported by LL::NG. It is a method to confirm a user's claimed identity by using a combination of two different factors between:

• something you know (login/password, . . . )

• something you have (U2F Key, TOTP, smartphone, mail, . . . )

• something you are (biometrics like fingerprints, . . . )

Applications Protection

• SAMLv2 protocol does not require a direct link between SP and IdP. All network exchanges go through the browser. • Handler code can be embedded by the protected applications but it requires to have access to the application code and to modify it. To avoid this, applications can be hidden behind Reverse Proxies that embed the Handler. So, protected applications are not exposed to users and directly requested (fig. 1).

Specific Features & Corresponding Plugins

AdaptativeAuthLevel A user reaches an authentication level depending on the selected authentication module, and eventually MFA. This plug-in adapts the authLevel depending on other conditions, like network, device, . . .

BruteForceProtection

Advanced Usages

• The ServiceToken (ST) handler is a mechanism implemented to protect the Server2Server exchanges. A web application (WebApp1) may need to request some other web applications (WebAppN) on behalf of the authenticated users. You just have to provide a header containing the ST to WebApp1. Then WebApp1 can request WebAppN by appending the ST in the 'X-LLNG-TOKEN' header that will be caught by the handler.

Presentation

  LemonLDAP::NG (LL::NG) is full AAA. It affords different ways for authenticating users, checks authorizations before accessing resources, and provides logs. (i) It relies on back-ends to store configuration and sessions. (ii) The Portal displays authentication screen, applications menu, and implements standard protocols. (iii) The Manager is the administration interface. (iv) The Handler protects applications working with HTTP headers.

Figure 1 :

 1 Figure 1: Main components

Figure 2 :

 2 Figure 2: Development cycle

Figure 3 :

 3 Figure 3: SAMLv2/Shibboleth kinematic

Figure 4 :

 4 Figure 4: OIDC Authorization Code Flow kinematic

Figure 5 :

 5 Figure 5: Handler kinematic

  To prevent brute force attack CheckDevOps Used for validating 'rules.json' files (fig. 7) CheckUser To check session attributes, access rights and transmitted headers for a specific user and URL ContextSwitching Admins can switch context other user CrowdSec A free and open-source security automation tool leveraging local IP behaviour detection and a communitypowered IP reputation system Impersonation Users can assume identity of another user Notification To notify some messages to users when log in REST/SOAP To provide corresponding services

Figure 6 :

 6 Figure 6: ServiceToken handler

Figure 7 :

 7 Figure 7: DevOps handler

Acknowledgements

Many thanks to ALEXANDRE KARIM from ESIEE Paris and the LemonLDAP::NG team for the figures and schemas.