Strongly Convex Optimization for Joint Fractal Feature Estimation and Texture Segmentation

Ho w Join t F ractal F eatures Estimation and T exture Segmen tation can b e cast in to a Strongly Con v ex Optimization Problem ? Barbara P ascal, Nelly Pustelnik and Patrice Abry ∗ † April 19, 2021 1 In tro duction Con text: unsup ervised texture segmentation. T extured images app ear as natur al mo dels for a large v ariety of real-w orld applications v ery different in nature. Often, fractal attributes are used to relev antly characterize such real- w orld textures. This is the case with biologic tissues [30], tomography-based pathology diagnostic [20, 29], art pain ting exp ertise [1], microfluidics [35], to name a few examples. Often in these applications, texture segmentation (i.e., partitioning images in to regions with homogeneous features) remains an on-going and ma jor chal- lenge. In computer vision or scene analysis, there exist numerous w ell-established and efficient metho ds to partition images mostly relying on their geometrical prop erties (cf. e.g., [19, 23, 31, 39]). F or textured images, segmentation remains c hallenging as geometry is more difficult to capture, relying essentially on the statistics of the texture features. Most classical texture segmen tation approaches rely on t wo-step procedures, consisting of first computing a priori chosen texture features and of, second, grouping features into regions with homogeneous features statistics. Attempts to impro ve these traditional approac hes address their to w main limitations: ar- bitrary a priori feature selection and sub optimalit y of the tw o-step pro cedure. F or instance, deep learning and neural netw ork approaches hav e strongly con- tributed to renew image segmen tation, av oiding notably the a priori selection of sp ecific features, and combining together feature estimation and segmentation (cf. e.g., [3]). Ho wev er, their use remains designed mostly for the specific con text of supervised segmentation: A (usually large) training dataset needs to be a v ail- able, together with large computer/memory p o wers. Often, suc h databases are not av ailable as expert annotations may be to o exp ensiv e (in time and money) to implement or techn ically b ey ond reac h. Besides tec hnical and financial is- sues, assembling such databases ma y allow raise non trivial ethical problems (cf. e.g., [13]). Along that line, a sev ere dra wback of neural netw ork approaches lies ∗ B. P ascal, N. Pustelnik and P . Abry are with the Universit ´ e de Ly on, ENS de Lyon, CNRS, Laboratoire de Physique, Lyon, F rance, e-mail: firstname.lastname@ens-lyon.fr † This work was supp orted b y Defi Imag’in SIR OCCO and by ANR-16-CE33-0020 Multi- F racs, F rance. 1 in the lack of interpr etability or even of identific ation of the features decisions are based on: Doctors for instance migh t legitimately remain reluctant to base diagnostic on non-iden tified features. Therefore, despite the p oten tial of deep learning, in contexts of absence of do cumen ted database, of imp ortance of accu- racy in estimated b oundaries, and of a need of understandabilit y , unsup ervised segmen tation strategies remain of critical importance. The present w ork is thus dedicated to the unsupervised segmentation of tex- tures in to regions well-c haracterized b y piecewise constan t fractal features (local regularit y and lo cal v ariance) formulated as a nonsmo oth strongly con v ex opti- mization problem, which can b e solv ed with fast and scalable algorithms. Related work. Amongst classical features commonly enrolled in texture seg- men tation, one can list e.g., Gab or or short-term F ourier transform co effi- cien ts [18, 22], fractal dimension [11], Amplitude/F requency Mo dulation mod- els [24]. More recently , fractal (or scale-free) features w ere also inv olved in tex- ture segmen tation (cf. e.g., [47]). Notably , lo c al r e gularity was used to quantify the fluctuations of roughness along the texture. Local regularity is quantified as an optimal lo cal p o wer-la w b ehavior across scales for some multiscale quan- tities [47]. Modulus of wa velet coefficients [27] were initially used as multiscale quan tities, follow ed by contin uous w av elet transform mo dulus maxima [28, 32]. More recen tly , wavelet le aders (lo cal suprema of discrete wa v elet co efficien ts), used here, were shown to permit a theoretically accurate and practically robust estimation of lo cal regularity [47], and successfully inv olv ed in texture charac- terization in e.g., [1, 29, 33]. W av elet leaders are used to estimate local regularity at each pixel mostly by linear regressions in log-log co ordinates. F or the segmentation step, local Gabor coefficient histograms were group ed using matrix factorization [48], ; textons were combined with brightness and color features to yield a multiscale contour detection pro cedure [4]. F urther, with fractal features, pixels sharing similar estimates of lo cal regularity are group ed together via a functional minimization strategy [38], using Rudin- Osher-F atemi (ROF) model [42] (i.e. total variation (TV) denoising). F urther, in [38], it was also attempted to combine b oth steps into a sin- gle one b y incorp orating the regression weigh t estimation in to the optimization pro cedure. The high computational burden implied b y iterativ ely solving opti- mization problems as well as by tuning the regularization parameter has ten- tativ ely been addressed by blo ck splitting appr o aches as suggested in [40] and explored in [34]. Strong conv exity [10] constitutes another recen tly prop osed trac k to address iterativ e optimization acceleration. Ho wev er, the prop osed joint features estimation and segmentation pro cedures, despite showing satisfactory and state-of-the-art segmentation performance suffered from ma jor limitations: Their high computing cost prev ents their use of large images and databases ; While based on a key fractal feature, lo cal regularity , they neglect c hanges in p o w er, a p oten tially relev ant information for texture segmen tation, notably to extract accurate region boundaries. Goals, contributions and outline. Aiming to address the abov e limitations, the ov erall goals of the present contribution are to inv estigate the p otential b enefits for texture segmentation brough t by • the use of joint estimates of fractal features (lo cal v ariance and lo cal reg- 2 ularit y) ; • the formulation of feature estimation and segmentation as a single step taking the form of a conv ex minimization ; • the deriv ation of fast and scalable iterativ e algorithms to solve the opti- mization problem, permitted by explicitly proving (and measuring) strong con vexit y of the data fidelity term in the ob jectiv e function ; • the comparisons of several optimization formulations, fa voring changes in features that are either co-lo calized in space or indep enden t ; • the deriv ation of an explicit sto c hastic pro cess, modeling piecewise fractal textures, th us p ermitting to conduct large-size Mon te-Carlo sim ulation for p erformance benchmarking and comparisons. T o that end, key concepts and definitions related to local regularit y , w av elet leaders and corresp onding state-of-the-art linear regression follow ed b y TV- based estimation procedures are recalled in Section 2. Tw o alternativ e ob jec- tiv e functions for the com bined estimation/segmen tation of lo cal v ariance and regularit y , resp ectively referred to as joint and c ouple d , are prop osed in Sec- tion 3, based on the same data fidelit y term but on tw o differen t total-v ariation based regularization strategies. Studying strong con vexit y and dualit y gaps, t wo classes of fast iterative algorithms ( Primal-Dual and Dual F orwar d-b ackwar d ), are devised. Piecewise homogeneous fractal textures are defined and studied in Section 4, th us enabling p erformance assessment from Mon te-Carlo sim ula- tions on well-con trolled syn thetic textures, as rep orted in Section 5. Conducting suc h simulations requires to address issues related to regularization parameter selection, iterativ e algorithm stopping criterion, and p erformance metrics. It p ermits relev an t answers to the final aim of assessing the actual b enefits of using b oth local v ariance and regularit y in texture segmentation with resp ect to issues such as sensitivity to fractal parameter c hanges, computational costs, impact of the differen t optimization form ulations. Finally , the p oten tial of the prop osed segmen tation approac hes is illustrated at w ork on real-w orld piecewise homogeneous textures taken from a publicly av ailable and documented texture database. Performance are also compared against a state-of-the-art segmenta- tion strategy . A ma tlab to olb o x implementing the analysis and syn thesis pro cedures de- vised here will b e made freely and publicly a v ailable at the time of publication. 2 Lo cal regularit y 2.1 Lo cal regularity and H¨ older exp onen t Lo cal regularity can classically b e assessed b y means of H¨ older exponent [21, 26], defined as follows: Definition 1 L et f : Ω → R denote a 2D r e al field define d on an op en set Ω ⊂ R 2 . The H¨ older exp onent h ( z 0 ) at lo c ation z 0 ∈ Ω is define d as the lar gest α > 0 such that ther e exists a c onstant χ > 0 , a p olynomial P z 0 of de gr e e lower 3 than α and a neighb orho o d V ( z 0 ) satisfying:  ∀ z ∈ V ( z 0 )    f ( z ) − P z 0 ( z )   ≤ χ k z − z 0 k α wher e k·k denotes the Euclidian norm. Definition 1 do es not provide a practical w ay to estimate lo cal regularity . Thus, the practical assessmen t of h ( z 0 ) usually relies on the use of multiscale quanti- ties, such as w av elet co efficien ts, or w av elet leaders [21, 26, 47]. 2.2 Lo cal regularity and w av elet leaders Because the practical aim is to analyse digitized images, definitions are giv en in a discrete setting. Let X = ( X n ) n ∈ Υ ∈ R | Υ | denote the digitized v ersion of the 2D real field f on a finite grid Υ = { 1 , . . . , N } 2 . Let d j = W j X denote the discrete w av elet transform (D WT) co efficien ts of X , at resolution j ∈ { j 1 , . . . , j 2 } , with W j : R N × N → R M j × M j the op erator form ulation of the DWT. F urther, let the wa velet leader L j,k , at scale 2 j and lo cation n = 2 j k , b e defined as the lo cal supremum of mo dulus of wa v elet co efficien ts in a small neigh b orho od across all finer scales, [21, 26, 47] L j,k = sup λ j 0 ,k 0 ⊂ 3 λ j,k | 2 j d j 0 ,k 0 | , with λ j,k = [ k 2 j , ( k + 1)2 j ) and 3 λ j,k = S p ∈{− 1 , 0 , 1 } 2 λ j,k + p . It w as prov en in [21, 26, 47] that w av elet leaders pro vide m ultiscale quan tities in trinsically tied to H¨ older exp onen ts, insofar as when X has H¨ older exponent h n at lo cation n , then the w av elet leaders L j,k satisfy: ( ∀ j ) L j,k ' η n 2 j h n as 2 j → 0 , with n = 2 j k ∈ λ j,k , (1) with η n prop ortional to lo cal v ariance of X around pixel n . 2.3 Lo cal regularity estimation 2.3.1 Linear regression Eq. (1) naturally leads to estimate h and v = log 2 η by means of linear regression in log-log co ordinates [21, 26, 47], denoted  b v LR , b h LR  : ( ∀ n ∈ Υ)  b v LR ,n b h LR ,n  = J − 1  S n T n  (2) with J =  R 0 R 1 R 1 R 2  and R m = X j j m , (3) and S n = X j log 2 L j,n , T n = X j j log 2 L j,n . (4) The sum P j implicitly stand for P j 2 j = j 1 with ( j 1 , j 2 ) the range of o ctav es in- v olved in the estimation. Though linear regressions pro vide estimates b oth for h and log 2 η , the later remains to date rarely used. Estimates of h w ere for instance used in [34, 38]. Ba yesian based extension of linear regressions were also and alternatively prop osed to estimate η and h [46]. 4 A sample of piecewise homogeneous fractal texture is displa yed in Fig. 1b, syn thesized using the mask sk etched in Fig. 1a. The corresponding linear re- gression based estimates b η LR ≡ 2 b v LR , and b h LR (Fig. 1(c-d)) sho w to o p o or p erformance (high v ariability) to p ermit to detect the tw o regions and the cor- resp onding boundary . (a) Mask (b) T exture X (c) b v LR (d) b h LR (e) b h ROF (f ) c M ROF Figure 1: Piecewise homogeneous fractal texture and lo cal v ariance and regularit y estimates. (a) Synthesis mask ( H 0 , Σ 0 ) = (0 . 5 , 0 . 6) and ( H 1 , Σ 1 ) = (0 . 8 , 0 . 65); (b) sample texture (see Section 4), (c) and (d) linear regression based estimates of lo cal v ariance and regularity b v LR and b h LR ; (e) and (f ) T otal v ariation based estimates of lo cal regularit y b h ROF and segmen ta- tion c M ROF obtained with Alg. 1. 2.3.2 T otal v ariation based estimates T o address the p oor estimation p erformance ac hieved b y linear regressions for lo cal estimation, it has b een proposed to denoise b h LR b y using an optimisation framew ork inv olving TV regularization, as a mean to enforce piecewise constant estimates [33, 38]: b h ROF = argmin h ∈ R | Υ | 1 2 k h − b h LR k 2 2 + λ TV ( h ) . (5) The abov e functional balances a data fidelit y term versus a total v ariation p e- nalization defined as: TV( h ) = k D h k 2 , 1 , (6) where the op erator D : R | Υ | → R 2 ×| Υ | , consists of horizon tal and vertical incre- men ts: ( D h ) n 1 ,n 2 =  h n 1 ,n 2 +1 − h n 1 ,n 2 h n 1 +1 ,n 2 − h n 1 ,n 2  (7) 5 and where the mixed norm, k·k 2 , 1 , defined for z = [ z 1 ; . . . ; z I ] > ∈ R I ×| Υ | as: k z k 2 , 1 = N 1 − 1 X n 1 =1 N 2 − 1 X n 2 =1 v u u t I X i =1 ( z i ) 2 n 1 ,n 2 (8) b oth ensures an isotropic and piecewise constan t behaviour of the estimates [38]. While showing substan tial b enefits compared to linear regression based es- timates and althought fast, the TV-based pro cedure ab o ve, prop osed in [38], suffers from several significan t limitations: i) It is restricted to the estimation of lo cal regularit y only and neglects local v ariance ; ii) It consists of a tw o-step pro cess (first apply linear regression to obtain b h LR , second apply TV to b h LR to obtain b h ROF ) and is hence potentially not optimal. Section 3 below will prop ose sev eral solutions that address these t wo limitations b y p erforming the estimate of b oth lo cal v ariance and regularity in a single pass at similar computational cost than the t wo-step procedure. 2.3.3 A p osteriori segmentation and lo cal regularit y global estima- tion When the num b er of regions is kno wn a priori , it has been sho wn in [7, 8] that a fast iterative thresholding p ost-processing procedure can b e applied to b h ROF to obtain a segmen tation of the TV denoised estimate. In [7, 8] were pro vided the- oretical guarantees linking this Threshold-ROF (T-R OF) strategy to Mumford- Shah-lik e segmentation. Alg. 1 explicitly describ es the pro cedure applied to b h ROF to obtain a t wo-region T-ROF segmen tation c M ROF with c M ROF , n = 0 if n ∈ Υ 0 and 1 otherwise. Figs. 1e, 1f rep ort examples of TV denoised estimate b h ROF and T-ROF segmen tation c M ROF . Algorithm 1: T-R OF: iterative thresholding of b h ROF Input: h Initialization: m [0] 0 = min n ∈ Υ h n , m [0] 1 = max n ∈ Υ h n . for t ∈ N ∗ do Compute the thr eshold: T [ t − 1] =  m [ t − 1] 0 + m [ t − 1] 1  / 2 Thr eshold h : Υ [ t ] 0 = { n | h n ≤ T [ t ] } , Υ [ t ] 1 = { n | h n > T [ t ] } Up date r e gion me an: m [ t ] 0 = 1 / | Υ 0 | X n ∈ Υ 0 h n , m [ t ] 1 = 1 / | Υ 1 | X n ∈ Υ 1 h n . Output: Υ 0 = Υ [ ∞ ] 0 , Υ 1 = Υ [ ∞ ] 1 F rom T-R OF segmentation, T exture X can b e interpreted as the concate- nation of tw o fractal textures X 0 and X 1 , each with its o wn uniform regularity and v ariance. Therefore, a p osteriori global estimates, b H 0 , ROF and b H 1 , ROF , for the lo cal regularity of eac h region, Υ 0 and Υ 1 , can b e obtained b y p erforming linear regressions applied to the logarithm of multiscale quan tities L j a veraged o ver Υ 0 (resp. Υ 1 ) (cf. [47]). 6 3 T otal v ariation based estimation of lo cal v ari- ance and regularit y 3.1 Design of the ob jectiv e function 3.1.1 Linear regression as functional minimization Linear regression can b e viewed as an optimization scheme. Indeed, setting v = log 2 η the following strictly conv ex functional (in v ariables ( v , h )): Φ( v , h ; L ) = 1 2 X j k v + j h − log 2 L j k 2 F ro , (9) has a unique minim um corresp onding to the linear regression estimates in 2. 3.1.2 P enalization T o improv e on the p o or p erformance of linear regression, a generic optimization framew ork can be proposed as ( b v , b h ) = arg min v ∈ R | Υ | , h ∈ R | Υ | Φ( v , h ; L ) + Ψ( v , h ) , where the data fidelity term Φ is that of the linear regression (cf.(9)), and where the p enalization term Ψ fav ors piecewise constancy of v and h . Two different strategies are prop osed: (i) Joint estimation : ( b v J , b h J ) = arg min v , h ∈ R | Υ | × R | Υ | Φ( v , h ; L ) + Ψ J ( v , h ) , (10) with Ψ J ( v , h ) = λ (TV( v ) + α TV ( h )) , (11) that do es not fav or changes in v and h that o ccur at same lo cation; (ii) Couple d estimation : ( b v C , b h C ) = arg min v , h ∈ R | Υ | × R | Υ | Φ( v , h ; L ) + Ψ C ( v , h ) , (12) with Ψ C ( v , h ) = λ TV α ( v , h ) , = λ    D v ; α D h  >   2 , 1 (13) where TV α couples spatial v ariations of v n and h n and thus fav or their o ccurrences at same lo cation. F or both constructions, λ > 0 and α > 0 constitute regularization hyperparam- eters that need to b e selected. F urther, the thresholding pro cedure in [7, 8] is generalized to “ joint ” and “ c ouple d ” estimation strategies: Alg. 1, applied to b h J (resp. b h C ) provides “ joint ” (resp. “ c ouple d ”) segmen tation c M J (resp. c M C ). Finally , from the tw o-region segmentation c M J (resp. c M C ) obtained with Alg. 1, one can use global techniques [47] to obtain a p osteriori estimates of the H¨ older exp onen ts of each region b H 0 , J and b H 1 , J (resp. b H 0 , C and b H 1 , C ) (cf. Sec. 2.3.2). 7 3.2 Minimization algorithms The data fidelity function Φ is Lipsc hitz differentiable. This is the case neither for Ψ J nor for Ψ C . Therefore, gradient descent metho ds are not appropriate to solv e (10) or (12). Instead, the use of proximal algorithms is adapted to the minimization of nonsmo oth functions [15]. While the non-differentiabilit y stems from the mixed norm k·k 2 , 1 , app earing b oth in Ψ J and in Ψ C , the linear op erator D in to the nonsmooth term makes the computation of the proximal op erator of total v ariation difficult [12, 37]. Below t wo algorithms are devised: the dual FIST A and the primal-dual. 3.2.1 Dual FIST A The w ell-known F ast Iterative Shrink age Thresholding Algorihtm (FIST A) al- gorithm [9] is here customized to Problems (10) and (12), to achiev e faster con vergences than with the basic dual forw ard-backw ard. Corresp onding itera- tions are detailed in Algorithms 2 and 3. Conv ergence guarantees are sp ecified in Theorems 1 and 2. Theorem 1 (Conv ergence of FIST A J ) The se quenc e ( v [ t ] , h [ t ] ) gener ate d by A lgorithm 2 c onver ges towar ds a solution of the joint estimation pr oblem (10) . Theorem 2 (Conv ergence of FIST A C ) The se quenc e ( v [ t ] , h [ t ] ) gener ate d by A lgorithm 3 c onver ges towar ds a solution of c ouple d estimation pr oblem (12) . Pro of 1 The pr o ofs of The or ems 1 and 2 stem dir e ctly fr om the choic e of de- sc ent p ar ameters pr op ose d in A lgs. 2 and 3 fol lowing the r e asoning in [9]. 3.2.2 Primal-dual Primal-dual algorithms [10, 16, 25, 44] can also b e customized to solve Prob- lems (10) and (12). This requires the deriv ation of a closed form expression for the proximal op erator asso ciated with the quadratic data fidelit y term Φ, provided in Propo- sition 1. Prop osition 1 (Computation of prox Φ ) F or every ( v , h ) ∈ R | Υ | × R | Υ | , de- noting ( p , q ) = prox Φ ( v , h ) ∈ R | Υ | × R | Υ | one has ( p = (1+ R 2 )( S + v ) − R 1 ( T + h ) (1+ R 0 )(1+ R 2 ) − R 2 1 , q = (1+ R 0 )( T + h ) − R 1 ( S + v ) (1+ R 0 )(1+ R 2 ) − R 2 1 . with S and T define d in (4) and R 0 , R 1 , R 2 define d in (3) . F urther, in [10], it was described ho w to tak e adv antage of the strong con- v exity 1 of the ob jective function to obtain fast implementation for primal-dual algorithms, with linear con vergence rates. In Proposition 2 b elo w, we pro ve that the data fidelit y term Φ is strongly con vex and deriv e the closed form expression of the strong con vexit y co efficien t. 1 A function ϕ : H → R , defined on an Hilbert space H , is said to be µ − strongly conv ex, for a given µ > 0 if the function y 7→ ϕ ( y ) − µ 2 k y k 2 2 is conv ex. When the function ϕ is differentiable, ϕ is µ − strongly con vex if and only if ( ∀ ( y , z ) ∈ H × H ) h∇ ϕ ( y ) − ∇ ϕ ( z ) , y − z i ≥ µ k y − z k 2 where h· , ·i denotes the Hilb ert scalar product and k·k the asso ciated scalar product. 8 Algorithm 2: FIST A J : Joint estimation (Pb. (10)) Initialization: Set u [0] ∈ R 2 ×| Υ | , ¯ u [0] = u [0] ; Set ` [0] ∈ R 2 ×| Υ | , ¯ ` [0] = ` [0] ; Let  S n , T n  defined in (4); Let J defined in (3); Set ( ∀ n )  v [0] n , h [0] n  > = J − 1  S n , T n  > ; Set b > 2 and τ 0 = 1; Set α > 0 and λ > 0; Set γ > 0 such that γ k J − 1 kk D k 2 < 1; for t ∈ N do Dual variable up date: u [ t +1] = pro x γ ( λ k . k 2 , 1 ) ∗  ¯ u [ t ] + γ D v [ t ]  ` [ t +1] = pro x γ ( λα k . k 2 , 1 ) ∗  ¯ ` [ t ] + γ D h [ t ]  FIST A p ar ameter up date τ t +1 = t + b b A uxiliary variable up date ¯ u [ t +1] = u [ t +1] + τ t − 1 τ t +1  u [ t +1] − u [ t ]  ¯ ` [ t +1] = ` [ t +1] + τ t − 1 τ t +1  ` [ t +1] − ` [ t ]  Primal variable up date  v [ t +1] h [ t +1]  =  v [ t ] h [ t ]  − J − 1  D ∗ ( u [ t +1] − u [ t ] ) D ∗ ( ` [ t +1] − ` [ t ] )  9 Algorithm 3: FIST A C : Coupled estimation (Pb. (12)) Initialization: Set u [0] ∈ R 2 ×| Υ | , ¯ u [0] = u [0] ; Set ` [0] ∈ R 2 ×| Υ | , ¯ ` [0] = ` [0] ; Let  S n , T n  defined in (4); Let J defined in (3); Set ( ∀ n )  v [0] n , h [0] n  > = J − 1  S n , T n  > ; Set b > 2 and τ 0 = 1; Set α > 0 and λ > 0; Set γ > 0 s. t. γ max(1 , α ) k J − 1 kk D k 2 < 1; for t ∈ N do Dual variable up date:  u [ t +1] ` [ t +1]  = pro x γ ( λ k . k 2 , 1 ) ∗ ¯ u [ t ] + γ D v [ t ] ¯ ` [ t ] + γ α D h [ t ] ! FIST A p ar ameter up date τ t +1 = t + b b A uxiliary variable up date ¯ u [ t +1] = u [ t +1] + τ t − 1 τ t +1  u [ t +1] − u [ t ]  ¯ ` [ t +1] = ` [ t +1] + τ t − 1 τ t +1  ` [ t +1] − ` [ t ]  Primal variable up date  v [ t +1] h [ t +1]  =  v [ t ] h [ t ]  − J − 1  D ∗ ( u [ t +1] − u [ t ] ) α D ∗ ( ` [ t +1] − ` [ t ] )  10 Prop osition 2 F unction Φ( v , h ; L ) is µ -str ongly c onvex w.r.t the variables ( v , h ) , with µ = χ wher e χ > 0 is the lowest eigenvalue of the symmetric and p ositive definite matrix J define d in Eq. (3) . Indeed, since ∇ Φ( v , h ; L ) = J ( v , h ) > is linear, the condition for Φ to b e µ − strongly conv ex can b e recasted as: h∇ Φ( v , h ; L ) , ( v , h ) i = D J ( v , h ) > , ( v , h ) > E ≥ µ k ( v , h ) k 2 (14) In tuitively a function with a large strong-con vexit y constan t µ is v ery stitched around its minim um, hence yielding a go o d theoretical conv ergence rate to ward the minimizer. The strong-conv exit y parameter is represented in Fig. 2 as a function of the range of scales inv olved in the estimation pro cedure (cf. (9)). Figure 2: Strong con vexit y constant as a function of the range of scales in volv ed in the estimation. The iterations of the devised fast primal-dual algorithms are detailed in Algorithms 4 and 5. Conv ergence guarantees are specified in Theorem 3 (resp. Theorem 4). Theorem 3 (Conv ergence of PD J ) The se quenc e ( v [ ` ] , h [ ` ] ) gener ate d by Al- gorithm 4 c onver ges towar ds a solution of the joint estimation pr oblem (10) . Theorem 4 (Conv ergence of PD C ) The se quenc e ( v [ ` ] , h [ ` ] ) gener ate d by Al- gorithm 5 c onver ges towar ds a solution of the joint estimation pr oblem (12) . Pro of 2 The pr o ofs of The or ems 3 and 4 c ome dir e ctly fr om the choic e of de- sc ent p ar ameters pr op ose d in A lgs. 4 and 5 fol lowing the r e asoning of [10]. 3.3 Dualit y gap T o ensure fair comparisons b et w een the differen t algorithms, w e construct a stopping criterion based on the duality gap . Let H , G Hilbert spaces, Θ : H → ] − ∞ , + ∞ ], Ξ : G → ] − ∞ , + ∞ ] and L : H → G a b ounded linear operator. F rom 11 Algorithm 4: PD J : Joint estimation (Pb. (10)) Initialization: Set v [0] ∈ R | Υ | , u [0] = D v [0] , ¯ u [0] = u [0] ; Set h [0] ∈ R | Υ | , ` [0] = D h [0] , ¯ ` [0] = ` [0] ; Set α > 0 and λ > 0; Set ( δ 0 , ν 0 ) such that δ 0 ν 0 k D k 2 < 1; for t ∈ N ∗ do Primal variable up date:  v [ t +1] h [ t +1]  = pro x δ t Φ  v [ t ] h [ t ]  − δ t D ∗ ¯ u [ t ] D ∗ ¯ ` [ t ] !! Dual variable up date: u [ t +1] = pro x ν t ( λ k·k 2 , 1 ) ∗  u [ t ] + ν t D v [ t +1]  ` [ t +1] = pro x ν t ( λα k·k 2 , 1 ) ∗  ` [ t ] + ν t D h [ t +1]  Desc ent steps up date: ϑ t = (1 + 2 µδ t ) − 1 / 2 , δ t +1 = ϑ t δ t , ν t +1 = ν t /ϑ t A uxiliary variable up date: ¯ u [ t +1] ¯ ` [ t +1] ! =  u [ t +1] ` [ t +1]  + ϑ t  u [ t +1] ` [ t +1]  −  u [ t ] ` [ t ]  Algorithm 5: PD C : Coupled estimation (Pb. (12)) Initialization: Set v [0] ∈ R | Υ | , u [0] = D v [0] , ¯ u [0] = u [0] ; Set h [0] ∈ R | Υ | , ` [0] = α D h [0] , ¯ ` [0] = ` [0] ; Set α > 0 and λ > 0 . Set ( δ 0 , ν 0 ) such that δ 0 ν 0 max(1 , α ) k D k 2 < 1; for t ∈ N ∗ do Primal variable up date:  v [ t +1] h [ t +1]  = pro x δ t Φ  v [ t ] h [ t ]  − δ t D ∗ ¯ u [ t ] α D ∗ ¯ ` [ t ] !! Dual variable up date:  u [ t +1] ` [ t +1]  = pro x ν t ( λ k . k 2 , 1 ) ∗  u [ t ] + ν t D v [ t +1] ` [ t ] + ν t α D h [ t +1]  Desc ent steps up date: ϑ t = (1 + 2 µδ t ) − 1 / 2 , δ t +1 = ϑ t δ t , ν t +1 = ν t /ϑ t A uxiliary variable up date: ¯ u [ t +1] ¯ ` [ t +1] ! =  u [ t +1] ` [ t +1]  + ϑ t  u [ t +1] ` [ t +1]  −  u [ t ] ` [ t ]  12 a primal optimization problem of the general form b x = arg min x ∈H Θ( x ) + Ξ( L x ), and its asso ciated dual problem b y = arg min y ∈G Θ ∗ ( − L ∗ y ) + Ξ ∗ ( y ), w e define: Γ( x ; y ) := Θ( x ) + Ξ( L x ) + Θ ∗ ( − L ∗ y ) + Ξ ∗ ( y ) , (15) Θ ∗ (resp. Ξ ∗ ), b eing the F enchel conjugate of Θ (resp. Ξ). Definition 2 The duality gap δ Γ asso ciate d to primal and dual optimization pr oblems is define d as the infimum: δ Γ := inf ( x , y ) ∈H×G Γ( x ; y ) . (16) Under some lo ose assumptions on Θ, Ξ and L , referred as str ong duality in [5], the duality gap δ Γ is zero and the resp ectiv e solutions b x and b y of primal and dual problems are c haracterized by Γ( b x ; b y ) = δ Γ = 0. Primal-dual algorithms presen ted ab o v e built a minimizing sequence ( x [ t ] , y [ t ] ) t ∈ N for Γ, therefore conv erging tow ards the unique p oin t achieving the infimum of Eq. (16): ( b x , b y ). Then the conv ergence of x [ t ] (resp. y [ t ] ) to wards b x (resp. b y ) can b e measured ev aluating the primal-dual functional Γ, as: Γ( x [ t ] ; y [ t ] ) − → t →∞ 0 . (17) The “ joint ” (10) and “ c ouple d ” (12) optimization problems, with primal v ariable x = ( v , h ) ∈ R | Υ | × R | Υ | , share the same data fidelity term Θ( x ) = Φ( v , h ; L ) and a p enalization based on the mixed norm Ξ = λ k·k 2 , 1 . Y et, the linear op erator L differ b etw een joint and c ouple d formulations, as: i) L x = [ D v , α D h ] ∈ R 2 × 2 | Υ | , for “ joint ” (10), ii) L x = [ D v ; α D h ] ∈ R 4 ×| Υ | , for “ c ouple d ” (12). F or Pb. (10) (resp. (12)), denoting y = ( u , ` ) the dual v ariable, the ev aluation of Γ J (resp . C) ( v , h ; u , ` ) requires computing: i) L ∗ (straigh tforward from D ∗ ), ii) Ξ ∗ = ( λ k·k 2 , 1 ) ∗ = ι k·k 2 , + ∞ ≤ λ (direct computation), iii) Θ ∗ = Φ ∗ whic h is devised in Prop osition 3. Prop osition 3 L et J and ( S , T ) define d in (3) and (4) , Φ ∗ ( v , h ; L ) = 1 2 h ( v , h ) > , J − 1 ( v , h ) > i + h ( S , T ) > , J − 1 ( v , h ) > i + C , (18) C = 1 2 h ( S , T ) > , J − 1 ( S , T ) > i − 1 2 X j (log 2 L j ) 2 . Pro of 3 By definition of the F enchel c onjugate, F ∗ ( v , h ; L ) = sup e v ∈ R | Υ | , e h ∈ R | Υ | h e v , v i + h e h , h i − F ( e v , e h ; L ) . (19) The supr emum is obtaine d at ( ¯ v , ¯ h ) such that, for every n ∈ Ω , ( v n − P j  ¯ v n + j ¯ h n − log 2 L j,n  = 0 h n − P j j  ¯ v n + j ¯ h n − log 2 L j,n  = 0 . (20) 13 or e quivalently,  R 0 ¯ v n + R 1 ¯ h n = v n + S n R 1 ¯ v n + R 2 ¯ h n = h n + T n (21) that yields to  ¯ v n ¯ h n  = J − 1  v n + S n h n + T n  (22) and it is then ne c essary to r e-inje ct this expr ession into the explicit expr ession of Ψ . 4 Piecewise homogeneous fractal textures Homogeneous fractal textures. Numerous mo dels of homogeneous textures w ere prop osed in the literature, mostly consisting of 2D extensions (referred to as fractional Brownian field fBf ) of fractional Brownian motion, the unique univ ariate Gaussian exactly selfsimilar pro cess with stationary incremen ts [14]. Ho wev er, in most applications, real world textures are b etter mo deled by sta- tionary pro cesses, hence b y increments of fBf, referred to as fractional Gaussian fields. F urther, to construct a formally relev an t definition for piecewise homo- geneous fractal textures, as originally prop osed here, it is easier to manipulate stationary pro cesses. F or Gaussian processes, space or frequency domain definitions of fBf are theoretically equiv alent and several v ariations w ere proposed in either domain, mostly consisting in different tuning of the 2D extension of the fractional in- tegration kernel underlying selfsimilar textures (cf. e.g., [2, 14, 36]). Y et, in practice, n umerical issues needs to b e accounted for, such as discrete sampling from a contin uous pro cess, or computation of in tegrals defined from a theoret- ically infinite support kernel and related b order effects. F urther, fast circulant em b edding matrix algorithms were prop osed, that yet do not p ermit to reach the full range 0 < H < 1 for the selfsimilarit y parameter [43]. Therefore, we will make use of a self-customized construction that combines an effective im- plemen tation scheme with excellen t theoretical and practical con trol of the local v ariance and regularity , while p ermitting direct extension to piecewise homoge- neous construction. F ollo wing e.g., [6, 17, 41], we start from an harmonizable representation of fBf B n , n ∈ R | Υ | , Υ = { 1 , . . . , N } 2 : B n = Σ C ( H ) Z e − i f .n − 1 k f k H +1 d e G ( f ) , (23) with d e G ( f ) the F ourier transform of a white Gaussian noise and C ( H ) = π 1 / 2 Γ( H + 1 / 2) 2 d/ 2 H Γ(2 H ) sin( π H )Γ( H + d/ 2) . The cov ariance function of the zero mean process B n reads: E  B n B m  = Σ 2 2  k n k 2 H + k m k 2 H − k n − m k 2 H  . (24) 14 Figure 3: Piecewise homogeneous fractal textures X generated using the mask display ed in Fig. 1a, with parameters ( H 0 , Σ 0 ) = (0 . 5 , 0 . 6) and ( H 1 , Σ 1 )) c hosen as: (a) (0 . 5 , 0 . 6), no c hange ; (b) (0 . 5 , 0 . 65), change in Σ ; (c) (0 . 8 , 0 . 6) c hange in H ; (d) (0 . 8 , 0 . 65), c hange in Σ and H . T exture Y n is constructed from incremen ts of B n as (with e 1 , e 2 , unitary v ectors in horizon tal and v ertical directions): Y n = Σ 2 δ H √ 1 − 2 H − 2 ( B n + δ e 1 − B n | {z } horizontal increment + B n + δ e 2 − B n | {z } vertical increment ) , (25) Prop osition 4 The field Y = ( Y n ) n ∈ Υ is a zer o-me an Gaussian pr o c ess with varianc e E h Y 2 n i = Σ 2 and c ovarianc e: E  Y n +∆ n Y n  = Σ 2 δ − 2 H 4 − 2 H ( k ∆ n + δ e 1 k 2 H + k ∆ n − δ e 1 k 2 H + k ∆ n + δ e 2 k 2 H + k ∆ n − δ e 2 k 2 H − 3 k ∆ n k 2 H − 1 2 k ∆ n + δ e 1 − δ e 2 k 2 H − 1 2 k ∆ n − δ e 1 + δ e 2 k 2 H ) . Piecewise homogeneous fractal textures. Piecewise homogeneous fractal textures X = ( X n ) n ∈ Υ are defined as the concatenation of M (pairwise disjoin t) regions, denoted (Υ m ) 0 ≤ m ≤ M − 1 , such that Υ = ∪ m Υ m , with textures in eac h region consisting of an homogeneous fractal Y ( m ) n (0 ≤ m ≤ M − 1), with v ariance and regularity (Σ 2 m , H m ): X n = Y ( m ) n , when n ∈ Υ m . Examples of piecewise fractal textures are shown in Fig. 3. 5 Estimation/Segmen tation p erformance 5.1 Mon te-Carlo sim ulation set-up The prop osed joint and c ouple d segmen tations are now compared in terms of p erformance and computational costs to the state-of-the-art T-ROF pro cedure, in the con text of a t wo-region segmentation, b y means of Monte-Carlo simula- tions. 5.1.1 Syn thetic textures Piecewise homogeneous fractal textures are generated as defined in Section 4, using the mask sho wn in Fig. 1a, with N = 512. They consist of t wo regions: i) a 15 bac kground, with v ariance and lo cal regularit y  Σ 2 0 , H 0  = (0 . 6 , 0 . 5) kept fixed, ii) a cen tral ellipse, for which v ariance and lo cal regularit y  Σ 2 1 , H 1  are v aried, as illustrated in Fig. 5. F or each configuration, 5 realizations of homogeneous fractal are generated and analyzed ; p erformance are rep orted as av erages ov er realizations. 5.1.2 W a v elet transform 2D-W a velet decompositions are p erformed using tensor pro duct wa velets con- structed from 1D-Daub ec hies orthonormal least asymmetric w av elets with N ψ = 3 v anishing momen ts [27]. W av elet leaders are computed as in Section 2.2. Lo- cal estimate b h LR in T-ROF is computed as in Sec. 2.3. Estimates from all pro cedures in volv e o cta ves ( j 1 , j 2 ) = (2 , 5). Octav e j 1 = 1 is a priori excluded as leaders are biased [45, 47]. 5.1.3 Optimization algorithm parameters T o achiev e b est conv ergence of the optimization schemes, descen t steps are c hosen as large as p ermitted: • Alg. 2 ( FIST A J ): γ = 0 . 99 /  k J − 1 kk D k 2  , • Alg. 3 ( FIST A C ): γ = 0 . 99 /  max(1 , α ) k J − 1 kk D k 2  , • Alg. 4 ( PD J ): δ 0 = ν 0 = 0 . 99 / k D k . • Alg. 5 ( PD C ): δ 0 = ν 0 = 0 . 99 / (max(1 , √ α ) k D k ), where k D k = 2 √ 2, with D defined in Eq. (7), and where k J − 1 k dep ends on the octav es in volv ed in estimation ( J defined in Eq. (3)): With ( j 1 , j 2 ) = (2 , 5), k J − 1 k ' 2 . 88. the inertia parameter is set to b = 4. 5.2 Issues in p erformance and algorithm comparison 5.2.1 Stopping criteria for proximal algorithms T o ensure fair comparisons either betw een algorithms or b etw een joint and c ou- ple d formulations (Pb. (10) vs. Pb. (12)), an effective stopping criterion is needed. Fig. 4 illustrates the conv ergence tow ard zero of tw o p oten tial can- didates, for the particular case, chosen as represen tative example, of the joint estimation using primal-dual Alg. 4, and for several c hoices of hyperparameters. First, the normalized increments of ob jectiv e function, as a function of the n umber of iterations, are observed (Fig. 4a) to decrease with an extremely ir- regular b eha vior, whic h makes them ill-suited to serv e as effective stopping criterion. Second, the normalize d primal-dual functional (with notations as in Sec. 3.3) e Γ [ t ] := Γ [ t ]   Θ( x [ t ] ) + Ξ( L x [ t ] )   +   Θ ∗ ( − L ∗ y [ t ] ) + Ξ ∗ ( y [ t ] )   , is observed (Fig. 4b) to decrease smoothly . Systematic insp ections of suc h de- creases together with that of the corresp onding achiev ed solutions lead us to 16 devise a stopping criterion, which is effectiv e b oth for primal-dual and forward- bac kward algorithms and for joint and c ouple d estimations, as, e Γ [ t ] • < 5 . 10 − 3 , for • ∈ { R OF , J } , and e Γ [ t ] C < 10 − 4 . (26) The use of the normalized quantities, e Γ [ t ] , makes this stopping criterion robust to v ariations of the h yp erparameters. F or practical purposes, w e also impose an upp er limit on the num b er of iterations: t < 2 . 5 10 5 . Figure 4: Stopping criteria for pro ximal algorithms. Con vergence of Alg. 4 solving joint problem (10) for fiv e pairs of hyperparameters ( λ, α ) ev al- uated with: (a) (Normalized) incremen ts of ob jectiv e function, (b) Normalized primal-dual functional e Γ [ t ] J . 5.2.2 Choice of regularization hyperparameters The choice of regularization parameters ( λ , α ) app earing in Pb (5), (10), (12) is of prime imp ortance as λ tunes the trade-off b et ween fidelity to the fractal mo del (1) and expected piecewise constancy , while α controls the relativ e w eigh t giv en to lo cal wa v elet log-v ariance v compared to lo cal regularity h , in the ( joint (11) and c ouple d (13)) total v ariation p enalization. The automated c hoice of the regularization parameters is a difficult issue, b ey ond the scope of the presen t work. In this study , a grid searc h strategy is used to find the parameters λ and α achieving the b est segmentation. In practice, logarithmically spaced ranges are used, from 10 − 1 to 10 3 for λ and from 10 − 2 to 10 3 for α . 5.2.3 P erformance assessmen t A natural p erformance criterion consists in comparing the achiev ed classifica- tion, denoted c M ROF , c M J and c M C resp ectiv ely for the three segmen tation pro cedures compared here, (ROF, joint and c ouple d ), to the mask in Fig 1a, regarded as ground truth. It leads to define the classific ation sc or e as the p er- cen tage of correctly lab eled pixels. Classification scores for c M ROF , c M J and 17 c M C , applied to different configurations of piecewise fractal textures are re- p orted in T able 1, together with the difference b et ween the a p osteriori global estimates obtained for eac h segmented regions Υ 0 and Υ 1 (cf. Section 2.3.3): d ∆ H • := b H 1 , • − b H 0 , • , for • ∈ { R OF , J , C } . 5.3 P erformance comparisons 5.3.1 Segmen tation and estimation p erformance Fig. 5 and T able 1 rep ort segmen tation and estimation performance for 7 differ- en t configurations and for the optimal set of hyperparameters (i.e., those that maximize the classification scores). Configurations I, I I I, V, VI corresp ond to a decrease in the difference b et w een the regularity of each region of the piecewise fractal texture: ∆ H = H 1 − H 0 , (hence to an increase in difficult y) for a fixed ∆Σ 2 = 0 . 1. While the segmenta- tion p erformance of the three pro cedures (T-R OF, T- joint and T- c ouple d ) are comparable for e asy configuration, those of T-ROF decrease drastically when ∆ H decreases while those of T- joint and T- c ouple d decrease significantly less. Along the same line, the estimation of ∆ H remains more s atisfactory at small ∆ H for T- joint and T- c ouple d than for T-ROF. It can also b e observed that the p erformance of T- c ouple d degrade slightly less than those T- joint . Configurations II, I II, IV corresp ond to a decrease in v ariance, ∆Σ 2 , (hence to an increase in difficult y) for a fixed ∆ H = 0 . 1 which can already b e regarded as a difficult case. As exp ected, T-R OF is not helped by the increase of v ari- ance b et w een IV and I I as estimation of lo cal regularity do es not dep end on v ariance [38, 47], and T-ROF segmen tation results are not satisfactory . It can also b e observ ed that performance of T- joint and T- c ouple d improv e when ∆Σ 2 increase, and again that the impro vemen t is sligh tly larger for T- c ouple d than for T- joint . These results p ermit to dra w tw o clear conclusions. First, there are quantifi- able b enefits in using the side information brought by ∆Σ 2 , notably when the c hanges in regularit y b ecome small (low ∆ H ): T- joint and T- c ouple d outp er- form T-ROF. Second, T- c ouple d , that, b y principle, fav or co-lo calized changes in regularity and v ariance, sho ws o verall better performance than T- joint , that do es not fav or co-lo calized changes. This is a satisfactory outcome as all the configurations c hosen follo w the a priori intuition, relev ant for real-world appli- cations, that changes of textures naturally imply co-lo calized c hanges in lo cal v ariance and lo cal regularit y . 5.3.2 Computational costs Comparisons in terms of computational costs b oth b et ween the three approac hes, and b etw een the tw o classes of proximal algorithms, dual forward-bac kw ard, standard and accelerated (FIST A), vs. primal-dual, standard and accelerated b y strong con vexit y (cf. Sec. 3), are reported in T able 2, for configurations I and I II (regarded as e asy and difficult ) considered as representativ e. Computational costs are rep orted in n umber of iterations actually used to reac h the stopping criterion and in real time, for the optimal set of hyperparameters and av eraged o ver 5 realizations. 18 (a) Explored ∆Σ 2 and ∆ H . I I I II I IV V VI ∆Σ 2 = 0 . 1 ∆ H = 0 . 2 ∆Σ 2 = 0 . 15 ∆ H = 0 . 1 ∆Σ 2 = 0 . 1 ∆ H = 0 . 1 ∆Σ 2 = 0 . 05 ∆ H = 0 . 1 ∆Σ 2 = 0 . 1 ∆ H = 0 . 05 ∆Σ 2 = 0 . 1 ∆ H = 0 . 025 T exture T-ROF T- joint T- c oupled Figure 5: Compared optimal segmentation. Piecewise fractal textures are c haracterized b y Σ 2 0 = 0 . 6, H 0 = 0 . 5 and different ∆Σ 2 , ∆ H as sketc hed in Fig. 1a. First ro w: T-R OF segmentation c M ROF . Second row: T- joint segmen- tation c M J . Third ro w: T- c ouple d segmentation c M C . I I I II I IV V VI ∆Σ 2 = 0 . 1 ∆ H = 0 . 2 ∆Σ 2 = 0 . 15 ∆ H = 0 . 1 ∆Σ 2 = 0 . 1 ∆ H = 0 . 1 ∆Σ 2 = 0 . 05 ∆ H = 0 . 1 ∆Σ 2 = 0 . 1 ∆ H = 0 . 05 ∆Σ 2 = 0 . 1 ∆ H = 0 . 025 T-ROF Score d ∆ H 86.7 ± 2.1% 0 . 21 ± 0 . 07 79.5 ± 1.2% 0 . 05 ± 0 . 02 78.5 ± 1.1% 0 . 05 ± 0 . 06 77.5 ± 2.9% 0 . 07 ± 0 . 04 69.9 ± 7.1% 0 . 01 ± 0 . 06 59.5 ± 2.4% 0 . 05 ± 0 . 07 T- joint Score d ∆ H 91.6 ± 1.7% 0 . 21 ± 0 . 06 91.5 ± 2.0% 0 . 07 ± 0 . 03 90.2 ± 1.9% 0 . 10 ± 0 . 02 84.2 ± 4.5% 0 . 04 ± 0 . 07 84.3 ± 3.2% 0 . 05 ± 0 . 02 74.7 ± 8.2% 0 . 11 ± 0 . 28 T- c ouple d Score d ∆ H 91.7 ± 1.7% 0 . 20 ± 0 . 05 91.9 ± 4.0% 0 . 06 ± 0 . 04 91.1 ± 1.5% 0 . 10 ± 0 . 02 85.5 ± 3.8% 0 . 08 ± 0 . 04 86.1 ± 4.3% 0 . 05 ± 0 . 02 74.3 ± 8.2% 0 . 06 ± 0 . 04 T able 1: Optimal segmentation p erformance for differen t configura- tions of fractal textures, av eraged ov er 5 realizations. Piecewise fractal textures are c haracterized by (Σ 2 0 , H 0 ) = (0 . 6 , 0 . 5) and differen t (∆Σ 2 , ∆ H ) as sk etched in Fig. 1a. First ro w: T-R OF segmentation. Second row: T- joint segmen tation. Third row: T- c ouple d segmentation. 19 T able 2 shows first that, as expected, accelerated algorithms are alw ays re- quire less iterations than non accelerated ones, thus generally leading to low er computational times (though this is not alwa ys the case with FIST A whose complexit y p er iteration is larger). Also, T-ROF sho ws alwa ys lo wer compu- tational costs compared to T- joint and T- c ouple d . This is expected as T-ROF only works with the regularity and do not use v ariance. FIST A vs. A c c eler ate d primal-dual. F or T-ROF, FIST A is o verall preferable to the accelerated primal-dual algorithm, as both show equiv alent computa- tional costs for Configuration I but FIST A is ten times faster (b oth in n umber of iterations and computation time) for Configuration I II. F or T- joint and T- c ouple d , for b oth configurations, accelerated primal-dual is faster than FIST A. F or T- c ouple d , in configuration I II, FIST A has actually not conv erged when meeting the the upp er limit of iterations. Therefore, FIST A is to b e preferred for T-ROF, while accelerated primal-dual algorithms are more relev an t for T- joint and T- c ouple d . T-joint vs. T-c ouple d. F ocusing on T- joint and T- c ouple d and th us on the accelerated primal-dual algorithm that is faster for these tw o methods, T able 2 sho ws that T- joint is solved 3 to 4 times faster (both in num b er of iterations and computational cost) than T- c ouple d . Configuration I Configuration I II T-ROF T- joint T- c ouple d T-ROF T- joint T- c ouple d Iterations (10 3 it.) DFB 96 ± 48 > 250 > 250 241 ± 18 > 250 > 250 FIST A 1 . 7 ± 0 . 4 50 . 2 ± 21 . 0 231 ± 37 3 . 7 ± 0 . 7 48 . 1 ± 3 . 4 > 250 PD 31 . 8 ± 17 . 0 > 250 > 250 201 ± 69 > 250 > 250 AcPD 1 . 5 ± 0 . 4 31 . 4 ± 4 . 6 125 ± 67 45 . 2 ± 43 40 . 5 ± 2 . 8 121 ± 42 Time (s) DFB 1 , 090 ± 520 4 , 840 ± 15 4 , 210 ± 76 2 , 010 ± 73 4 , 810 ± 215 4 , 200 ± 76 FIST A 16 ± 4 1 , 030 ± 410 4 , 800 ± 560 30 ± 5 989 ± 64 5 , 110 ± 340 PD 297 ± 150 4 , 180 ± 69 4 , 110 ± 43 1 , 580 ± 490 4 , 150 ± 18 4 , 100 ± 15 AcPD 15 ± 4 619 ± 96 2 , 420 ± 1 , 300 349 ± 330 785 ± 59 2 , 320 ± 790 T able 2: Number of iterations and computational time necessary to reach Condition (26) for the different proximal algorithms inv estigated, illustrated on tw o configurations I (∆ H = 0 . 2, ∆Σ 2 = 0 . 1) and I II (∆ H = 0 . 1, ∆Σ 2 = 0 . 1). DFB : Dual F orward-Bac kw ard, FIST A : inertial acceleration of DFB, PD : primal-dual, AcPD : strong-con vexit y based acceleration of PD. 5.3.3 Ov erall comparison As an o verall conclusion, results rep orted ab ov e sho w that there are b enefits to use together lo cal regularit y and v ariance, compared to regularit y only , when the c hanges in regularity become small. This implies switching from accelerated dual forw ard-backw ard algorithms (FIST A, for T-ROF) to accelerated primal- dual algorithms (for T- joint and T- c ouple d ). F or difficult configurations, T- c ouple d (sligh tly) outp erforms T- joint in terms of segmen tation p erformance, at the price of non negligible increases of compu- tational costs. In that sense, T- joint can be considered a reasonable trade-off b et w een too p oor segmentation p erformance (as those of T-R OF) and to o large computational costs (as those of T- c ouple d ). 20 5.4 Real-w orld textures The actual p oten tial of the prop osed segmentation strategies is further illus- trated on real-w orld textures, arbitrarily c hosen from the commonly used large, do cumen ted and publicly av ailable University of Maryland, HighR esolution (UMD HR) 2 texture dataset. This dataset consists of 50 classes of homogeneous tex- tures, eac h consisting of 40 images of same textures under differen t imaging conditions (ligh t, angle,. . . ). F rom arbitrarily chosen pairs of images, piecewise homogeneous textures are constructed b y including textures one into the other using the ellipse-shaped mask used ab o v e, as sketc hed and illustrated in Fig. 6. Homogeneous textures are centered and normalized in v ariance indep enden tly and prior to inclusion, to av oid that detection could be based on sole c hange of means or v ariances. The T-R OF, T- Joint and T- Couple d segmentation approaches prop osed here are applied to and compared on these real-world piecewise homogeneous tex- tures. F urther, they are compared against one of the state-of-the-art texture segmen tation proc edure, based on matrix factorization and lo cal sp ectral his- tograms clustering, proposed in [48] and b eing a joint feature selection and segmen tation procedure, and hereafter referred to as LSHC-MF . Fig. 7 rep orts and compares segmentation p erformance, in terms of p er- cen tage of pixels correctly lab elled. F or the T-ROF, T- Joint and T- Couple d approac hes, optimal p erformance are rep orted, after a grid search for the se- lection of the hyperparameters λ and α . F or LSHC-MF , optimal p erformance are also rep orted, after a grid searc h on the window size, within whic h local histogram statistics are computed. Fig. 7 clearly sho ws that the LSHC-MF approach fails to segment correctly the piecewise homogeneous texture and is outp erformed by the T-ROF, T- Joint and T- Couple d strategies. It also shows that T-ROF, while not b eing inconsistan t, does not yield satisfactory segmen tation, whereas T- Joint and T- Couple d do. F urther, it is worth noting that not only the Couple d strategy reac hes (sligh tly) b etter segmentation performance compared to the Joint one, but that it also yields smo other b oundaries, an outcome of utmost relev ance for decision/diagnostic in numerous real-world applications, cf. e.g. [35]. Suc h significan t improv emen t in segmentation p erformance is at the cost tough of significan tly higher computational cost: When LSHC-MF requires of the order of 1s to p erform its most accurate segmen tation, as in Fig. 7(d), The T-R OF, T- Joint and T- Couple d segmentation approac hes require of the orders of 10s, 600s 1800s respectively , in agreement with computational cost ev aluation rep orted in T ab. 2. Consisten t conclusions can b e dra wn for other c hoices of classes of textures in the UMD HR dataset. 6 Conclusion and future work The present article has significantly adv anced the state-of-the-art in the seg- men tation of piecewise fractal textures. First, it has b een proposed to base the segmen tation of fractal textures not only on the estimation of the sole lo cal regularity parameter, but to use an additional lo cal parameter, the log-wa v elet v ariance, tightly related to the lo cal 2 http://legacydirs.umiacs.umd.edu/ ∼ fer/website-texture/texture.htm 21 (a) T exture A (b) T exture B (c) Mask (d) Mixture A - B Figure 6: Real world textures. t wo samples of real-world textures taken from the UMD HR texture dataset ((a) and (b)), piecewise constant mask (c) and piecewise homogeneous texture (d). 82 . 6% 86 . 1% 87 . 9% 55 . 2% (a) T-ROF (b) Joint (c) Coupled (d) LSHC-MF Figure 7: Segmentation p erformance on real-world textures. F or the piecewise homogeneous texture sho wn in Fig. 6(d), p erformance are computed in term of the p ercentage of well-classified pixels for T-ROF (a), Joint (b), Couple d (c) and LSHC-MF (d). 22 v ariance of the textures. Tw o v ariations were in vestigated, c ouple d and joint , that respectively enforce or not co-lo calized c hanges in regularit y and v ariance. It has b een shown, using large size Monte Carlo simulations, that the use of this additional features impro ves drastically segmentation performance when the difference in regularit y b ecomes negligible. This yet comes at the price of a non negligible increase in computational costs. Therefore, a second contribution has b een to construct accelerated primal- dual algorithms, requiring the explicit calculation of the strong conv exity con- stan t underlying the data fidelit y term form. The ac hiev ed substan tial reduction in computational costs has turned critical both to b e able to conduct large size Mon te Carlo sim ulation and to p erform the greedy search of the optimal set of h yp er–parameters. This low computational cost is also crucial for application on real-world data. The inv estigations reported here hav e permitted to show that accelerated primal-dual algorithms outp erform accelerated dual forward-bac kward (FIST A- t yp e) algorithms for piece wise fractal texture segmentation as so on as the joint use of regularit y and v ariance is required. F urther, they sho wed that the c ouple d form ulation, that fa vor co-lo calized c hanges in regularity and v ariance, p erforms b etter than the joint formulation, yet at the price of a significan tly larger com- putational cost. Thus, depending on budget constraints on time and requested qualit y of the solution, the joint formulation can b e regarded as an effective trade-off. The prop osed theoretical form ulations for piecewise fractal texture segmen- tation and the corresp onding accelerated algorithms are matured enough for applications on real-w orld data and competitive with resp ect to state-of-the-art as shown in Sec. 5.4. Application to the segmen tation of m ultiphasic flo ws is under curren t inv estigations. The automation of the tuning of the h yp erpa- rameters, whic h is a crucial p oin t for real-world applications, as underlined in Sec. 5.4, is also being in vestigated. Extensions to piecewise m ultifractal textures are further targeted. A ma tlab to olb o x implementing the analysis and syn thesis pro cedures de- vised here will b e made freely and publicly a v ailable at the time of publication. References [1] P . Abry , S. Jaffard, and H. W endt. When Van Gogh meets Mandelbrot: Multifractal classification of painting’s texture. Signal P r o c ess. , 93(3):554– 572, 2013. [2] P . Abry and F. Sellan. The w av elet-based syn thesis for fractional Brownian motion prop osed by F. Sellan and Y. Mey er: Remarks and fast implemen- tation, 1996. [3] V. Andrearczyk and P . F. Whelan. Using filter banks in conv olutional neural netw orks for texture classification. Pattern R e c o gn. L ett. , 84:63–69, 2016. [4] P . Arb elaez, M. Maire, C. F o wlkes, and J. Malik. Contour detection and hierarc hical image segmen tation. IEEE T r ans. Pattern Anal. Match. Int. , 33(5):898–916, 2011. 23 [5] H. H. Bauschk e and P . L. Com b ettes. Convex Analysis and Monotone Op er ator The ory in Hilb ert Sp ac es . Springer, New Y ork, 2011. [6] H. Bierm´ e, M. Meerschaert, and H. Scheffle. Op erator scaling stable ran- dom fields. Sto chastic Pr o c esses and their Applic ations , 117(3):312–332, 2007. [7] X. Cai, R. Chan, C.-B. Schonlieb, and T. Steidl, G.and Zeng. Link age b et w een Piecewise Constant Mumford-Shah mo del and ROF model and its virtue in image segmen tation. arXiv pr eprint arXiv:1807.10194 , 2018. [8] X. Cai and G. Steidl. Multiclass segmentation by iterated ROF thresh- olding. In International Workshop on Ener gy Minimization Metho ds in Computer Vision and Pattern R e c o gnition , pages 237–250. Springer, 2013. [9] A. Cham b olle and C. Dossal. On the conv ergence of the iterates of “FIST A”. J. Optim. The ory Appl. , 166(3):25, August 2015. [10] A. Chambolle and T. Pock. A first-order primal-dual algorithm for conv ex problems with applications to imaging. J. Math. Imag. Vis. , 40(1):120–145, 2011. [11] B. B. Chaudhuri and N. Sark ar. T exture segmentation using fractal dimen- sion. IEEE T r ans. Pattern Anal. Match. Int. , 17(1):72–77, 1995. [12] C. Chaux, P .L. Com b ettes, and V.R. Pesquet, J.-C. W a js. A v ariational for- m ulation for frame-based inv erse problems. Inverse Pr oblems , 23(4):1495– 1518, Jun. 2007. [13] M. Cimp oi, S. Ma ji, I. Kokkinos, and A. V edaldi. Deep filter banks for texture recognition, description, and segmen tation. Int. J. Comp. Vis. , 118(1):65–94, 2016. [14] S. Cohen and J. Istas. F r actional fields and applic ations . Springer, 2013. [15] P . L. Combettes and J.-C. Pesquet. Pro ximal splitting metho ds in signal pro cessing. In H. H. Bauschk e, R. S. Burac hik, P . L. Com b ettes, V. Elser, D. R. Luke, and H. W olk owicz, editors, Fixe d-Point Algorithms for Inverse Pr oblems in Scienc e and Engine ering , pages 185–212. Springer-V erlag, New Y ork, 2011. [16] L. Condat. A primal-dual splitting metho d for conv ex optimization in volv- ing lipsc hitzian, proximable and linear composite terms. J. Optim. The ory Appl. , 158(2):460–479, 2013. [17] G. Didier, M. M. Meerschaert, and V. Pipiras. Domain and range symme- tries of op erator fractional Brownian fields. Sto chastic Pr o c esses and their Applic ations , 128(1):39–78, 2018. [18] D. Dunn, W. E. Higgins, and J. W akeley . T exture segmentation using 2- D Gabor elementary functions. IEEE T r ans. Pattern Anal. Match. Int. , 16(2):130–149, 1994. [19] P . Getreuer. Chan-Vese segmen tation. Image Pr o c essing On Line , 2:214– 224, 2012. 24 [20] M. A. Ibrahim, O. A. Ojo, and P . A. Olu wafiso ye. Identification of emphy- sema patterns in high resolution computed tomography images. Journal of Biome dic al Engine ering and Informatics , 4(1):16, 2017. [21] S. Jaffard. W a velet techniques in multifractal analysis. F r actal Ge ometry and Applic ations: A Jubile e of Beno ˆ ıt Mandelbr ot, M. L apidus and M. van F r ankenhuijsen Eds., Pr o c e e dings of Symp osia in Pur e Mathematics , 72(2):91–152, 2004. [22] A. K. Jain and F. F arrokhnia. Unsupervised texture segmen tation using Gab or filters. Pattern R e c o gn. , 24(12):1167–1186, 1991. [23] M. Jung, G. P eyr´ e, and L. D. Cohen. Non-lo cal Activ e Con tours. In Pr o c. SSVM , volume 6667, pages 255–266, Ein-Gedi, Israel, 2011. Springer. [24] I. Kokkinos, G. Ev angelop oulos, and P . Maragos. T exture analysis and segmen tation using mo dulation features, generative mo dels, and weigh ted curv e evolution. IEEE T r ans. Pattern Anal. Match. Int. , 31(1):142–157, 2009. [25] N. Komodakis and J.-C. Pesquet. Playing with duality: An ov erview of re- cen t primal-dual approaches for solving large-scale optimization problems. IEEE Signal Pr o c ess. Mag. , 32(6):31–54, No v. 2015. [26] R. Leonarduzzi, H. W endt, P . Abry , S. Jaffard, C. Melot, S. G. Roux, and M. E. T orres. p-exp onen t and p-leaders, Part I I: Multifractal Analysis. Re- lations to Detrended Fluctuation Analysis. Physic a A , 448:319–339, 2016. [27] S. Mallat. A wavelet tour of signal pr o c essing . Academic Press, San Diego, USA, 1997. [28] S. Mallat and S. Zhong. Characterization of signals from multiscale edges. IEEE T r ans. Pattern A nal. Match. Int. , 14(7):710–732, 1992. [29] Z. Marin, K. A. Batchelder, B. C. T oner, L. Guimond, E. Gerasimov a- Chec hkina, A. R. Harrow, A. Arneodo, and A. Khalil. Mammographic ev- idence of micro en vironmen t c hanges in tumorous breasts. Me dic al Physics , 44(4):1324–1336, 2017. [30] M. T. McCann, D. G. Mixon, M. C. Fickus, C. A. Castro, J. A. Ozolek, and J. Kov acevic. Images as o cclusions of textures: A framework for seg- men tation. IEEE T r ans. Image Pr o c ess. , 23(5):2033–2046, 2014. [31] D. Mumford and J. Shah. Optimal appro ximations b y piecewise smo oth functions and asso ciated v ariational problems. Comm. Pur e Applie d Math. , 42(5):577–685, 1989. [32] J.F. Muzy , E. Bacry , and A. Arneo do. W av elets and m ultifractal formalism for singular signals: Application to turbulence data. Phys. R ev. L ett. , 67(25):3515–3518, 1991. [33] J.D.B. Nelson, C. Nafornita, and A. Isar. Semi-local scaling exp onen t estimation with b o x-p enalt y constraints and total-v ariation regularization. IEEE T r ans. Image Pr o c ess. , 25(7):3167–3181, 2016. 25 [34] B. Pascal, N. Pustelnik, P . Abry , and J.-C. Pesquet. Blo c k-co ordinate pro ximal algorithms for scale-free texture segmen tation. In Pr o c. Int. Conf. A c oust., Sp e e ch Signal Pr o c ess. , Calgary , Alb erta, Canada, Apr. 15-20 2018. [35] B. Pascal, N. Pustelnik, P . Abry , M. Serres, and V. Vidal. Join t Estimation of Lo cal V ariance and Lo cal Regularity for T exture Segmentation. Appli- cation to Multiphase Flow Characterization. In Pr o c. Int. Conf. Image Pr o c ess. , pages 2092–2096, Athens, Greece, 2018. IEEE. [36] R. M. P ereira, C. Garban, and L. Chevillard. A dissipativ e random v elo c- it y field for fully dev elop ed fluid turbulence. Journal of Fluid Me chanics , 794:369–408, 2016. [37] N. Pustelnik, C. Chaux, and J.-C. P esquet. P arallel ProXimal Algorithm for image restoration using h ybrid regularization. IEEE T r ans. Image Pr o- c ess. , 20(9):2450–2462, 2011. [38] N. Pustelnik, H. W endt, P . Abry , and N. Dobigeon. Com bining lo cal reg- ularit y estimation and total v ariation optimization for scale-free texture segmen tation. IEEE T r ans. Computational Imaging , 2(4):468–479, 2016. [39] S. Ray and R. H. T uri. Determination of num b er of clusters in K-means clustering and application in colour image segmen tation. In Pr o c. Int. Conf. A dv. Pattern R e c o gnition and Digital T e chniques , pages 137–143. Calcutta, India, 1999. [40] A. Rep etti, E. Chouzenoux, and J.-C. Pesquet. A parallel blo ck-coordinate approac h for primal-dual splitting with arbitrary random blo c k selection. In Pr o c. Eur. Sig. Pr o c. Confer enc e , pages 235–239, Nice, F rance, Aug. 31-Sept. 4, 2015. [41] S.G. Roux, M. Clausel, B. V edel, S. Jaffard, and P . Abry . Self-Similar Anisotropic T exture Analysis: the Hyp erbolic Wav elet transform contribu- tion. IEEE T r ans. Image Pr o c ess. , 22(11):4353 – 4363, 2013. [42] L. I. Rudin, S. Osher, and E. F atemi. Nonlinear total v ariation based noise remo v al algorithms. Physic a D: Nonline ar Phenomena , 60(1-4):259–268, 1992. [43] M. L. Stein. F ast and exact simulation of fractional Brownian surfaces. Journal of Computational and Gr aphic al Statistics , 11(3):587–599, 2002. [44] B.C. V ˜ u. A splitting algorithm for dual monotone inclusions inv olving co coercive operators. A dv. Comput. Math. , 38(3):667–681, Apr. 2013. [45] H. W endt, P . Abry , and S. Jaffard. Bo otstrap for Empirical Multifractal Analysis. IEEE Signal Pr o c ess. Mag. , 24(4):38–48, 2007. [46] H. W endt, S. Combrexelle, Y. Altmann, J.-Y. T ourneret, S. McLaughlin, and P .Abry . Multifractal analysis of m ultiv ariate images using gamma Mark ov random field priors. SIAM J. Imaging Sci. , 11(2):1294–1316, 2018. [47] H. W endt, S. G. Roux, P . Abry , and S. Jaffard. W av elet leaders and bo ot- strap for m ultifractal analysis of image s. Signal Pr o c ess. , 89(6):1100–1114, 2009. 26 [48] J. Y uan, D. W ang, and A. M. Cheriyadat. F actorization-based texture segmen tation. IEEE T r ans. Image Pr o c ess. , 24(11):3488–3497, Nov em b er 2015. 27