Year: 2021

Abstract:
State-of-the-art maximum entropy models for texture synthesis are built from statistics relying on image representations defined by convolutional neural networks (CNN). Such representations capture rich structures in texture images, outperforming wavelet-based representations in this regard. However, conversely to neural networks, wavelets offer meaningful representations, as they are known to detect structures at multiple scales (e.g. edges) in images. In this work, we propose a family of statistics built upon non-linear wavelet based representations, that can be viewed as a particular instance of a one-layer CNN, using a generalized rectifier non-linearity. These statistics significantly improve the visual quality of previous classical wavelet-based models, and allow one to produce syntheses of similar quality to state-of-the-art models, on both gray-scale and color textures. We further provide insights on memorization effects in these models.

heading: 1 INTRODUCTION
text: Textures ares spatially homogeneous images, consisting of similar patterns forming a coherent ensemble. In texture modeling, one of the standard approaches to synthesize textures relies on defining a maximum entropy model (Jaynes, 1957) using a single observed image (Raad et al., 2018). It consists of computing a set of prescribed statistics from the observed texture image, and then generating synthetic textures producing the same statistics as the observation. If the statistics correctly describe the structures present in the observation, then any new image with the same statistics should appear similar to the observation. A major challenge of such methods resides in finding a suitable set of statistics, that can generate both high-quality and diverse synthetic samples. This problem is fundamental as it is at the heart of many texture related problems. For example, in patch re-arrangement methods for texture modeling, these statistics are used to compute high-level similarities of image patches (Li & Wand, 2016; Raad et al., 2018). Such statistics are also used in texture interpolation for probing visual perception (Vacher et al., 2020), style transfer and image inpainting (Gatys et al., 2016; Laube et al., 2018).
A key question along this line of research is to find what it takes to generate natural textures. This problem was originally posed in Julesz (1962), in which the author looks for a statistical characterization of textures. In the classical work of Portilla & Simoncelli (2000), the authors presented a model whose statistics are built on the wavelet transform of an input texture image. These statistics were carefully chosen, by showing that each of them captured a specific aspect of the structure of the image. This model produces satisfying results for a wide range of textures, but fails to reproduce complex geometric structures present in some natural texture images. Figure 1 presents a typical example composed of radishes, and synthetic images from three state-of-the-art models developed over the last few decades. To address this problem, the work of Gatys et al. (2015) proposes to use statistics built on the correlations between the feature maps of a deep CNN, pre-trained on the ImageNet classification problem (Deng et al., 2009; Simonyan & Zisserman, 2014). While this model produces visually appealing images, these statistics are hard to interpret. The work of Ustyuzhaninov et al. (2017) made a significant simplification of such statistics, by using the feature maps of a one-layer rectifier CNN with random filters (without learning). A crucial aspect of this simplification relies on using multi-scale filters, which are naturally connected to the wavelet transform. In this paper, we propose a wavelet-based model, more interpretable than CNN-based models (with learned or random filters), to synthesize textures with complex geometric structures. It allows to bridge the gap between the classical work of Portilla & Simoncelli (2000), and state-of-the-art models.
This model is built on the recent development of the phase harmonics for image representations and non-Gaussian stationary process modeling (Mallat et al., 2020; Zhang & Mallat, 2021). The phase harmonics are non-linear transformations that adjust the phase of a complex number. In Portilla & Simoncelli (2000); Zhang & Mallat (2021), the authors illustrate that the phase dependencies between wavelet coefficients across scales contain important information about the geometric structures in textures and turbulent flows, and that they can be captured by applying the phase harmonics to complex wavelet coefficients. Remarkably, Mallat et al. (2020) show that the phase harmonics admit a dual representation, closely related to the rectifier non-linearity in CNNs.
Our main contributions are:
• We develop a family of texture models based on the wavelet transform and a generalized rectifier non-linearity, that significantly improves the visual quality of the classical waveletbased model of Portilla & Simoncelli (2000) on a wide range of textures. It relies on introducing spatial shift statistics across scales to capture geometric structures in textures.
• By changing the number of statistics in our models, we show explicitly the trade-off on the quality and diversity of the synthesis. When there are too many statistics, our model tends to memorize image patches. We further investigate such memorization effects on non-stationary images and find that it sometimes relies on what statistics are chosen, rather than on how many.
• Through the modeling of geometric structures in gray-scale textures, our model indicates the possibility of reducing significantly the number of statistics in the works of Gatys et al. (2015) and Ustyuzhaninov et al. (2017), to achieve a similar visual quality.
The rest of the paper is organized as follows: Section 2 reviews the framework of microcanonical maximum-entropy models, build upon a general family of covariance statistics. We then present our model for both gray-scale and color textures in Section 3. Section 4 shows synthesis results of our model, compared with state-of-the-art models. Finally, in Section 5, we discuss possible improvements of our model.
Notations Throughout the paper,N denotes a positive integer. A gray-scale image x is an element of RN×N , i.e. x = x(u), u ∈ ΩN , with ΩN := {0, · · · , N−1}2. A color image x = {xc}c=1,2,3 is an element of R3×N×N , or equivalently, each xc ∈ RN×N . We shall denote x̄ the observed texture (observation), which is assumed to be a realisation a random vector X . For any complex number z ∈ C, z∗ is the complex conjugate of z, Real(z) its real part, |z| its modulus, and ϕ(z) its phase.

heading: 2 MICROCANONICAL COVARIANCE MODELS
text: We briefly review the standard framework of micro-canonical maximum-entropy models for textures. To reliably estimate the statistics in these models, we assume that a texture is a realization of a stationary and ergodic process X (restricted to ΩN ). We then review a special family of statistics that are used in the state-of-the-art texture models (mentioned in Figure 1), based on covariance statistics of an image representation.
1As the model from Ustyuzhaninov et al. (2017) uses on random filters, we shall use the abbreviation RF.

heading: 2.1 FRAMEWORK
text: Given a observation texture x̄, we aim at generating new texture images, similar but different from x̄. To that end, a classical method is to define a set of statistics Cx̄, computed on the observation, and try to sample from the microcanonical set
{x : ‖Cx− Cx̄‖ ≤ }, where ‖ · ‖ denotes the L2 norm. Under the stationary and ergodic assumption of X , one can construct Cx as a statistical estimator of E(CX), from a complex-valued representation Rx.2 The set of covariance statistics Cx of a model can then be constructed by computing an averaging over the spatial variable u, i.e.
Cx(γ, γ′, τ) := 1 |ΩN | ∑ u∈ΩN Rx(γ, u)Rx(γ′, u− τ)∗, (1)
for (γ, γ′, τ) ∈ Υ ⊆ Γ × Γ × ΩN . The statistics Cx(γ, γ′, τ) can be interpreted as estimating the covariance (resp. correlations) betweenRX(γ, u) andRX(γ′, u−τ) for zero-meanRX (resp. nonzero meanRX). The ergodicity assumption ensures that whenN is large enough, the approximation Cx̄ ' E(CX) over Υ should hold with high probability. Under these conditions, it makes sense to sample the microcanonical set in order to generate new texture samples.
This framework encompasses a wide range of state-of-the-art texture models3. In particular, the PS model developed in Portilla & Simoncelli (2000) takes inspiration from the human early visual system to define a multi-scale representation based on the wavelet transform of the image (Heeger & Bergen, 1995). We next review a family of covariance model which generalizes the statistics in the PS model. We write CM the statistics for a specific model M that uses the representation RM.

heading: 2.2 WAVELET PHASE HARMONIC COVARIANCE MODELS
text: We review a family of microcanonical covariance models defined by a representation built upon the wavelet transform and phase harmonics. It defines a class of covariance statistics that capture dependencies between wavelet coefficients across scales.

heading: 2.2.1 WAVELET TRANSFORM
text: The wavelet transform is a powerful tool in image processing to analyze signal structures, by defining a sparse representation (Mallat, 2001). For texture modeling, we consider oriented wavelets to model geometric structures in images at multiple scales. They include the Morlet wavelets and steerable wavelets, proposed in Goupillaud et al. (1984); Simoncelli & Freeman (1995); Unser & Chenouard (2013). In particular, the Simoncelli steerable wavelets have been used to model a diverse variety of textures in Portilla & Simoncelli (2000).
Oriented wavelets are defined by the dilation and rotation of a complex function ψ : R2 7→ C on a plane. Let rθ denote the rotation by angle θ in R2. They are derived from ψ with dilations by factors 2j , for j ∈ {0, 1, · · · , J − 1}, and rotations rθ over angles θ = `π/L for 0 ≤ ` < L, where L is the number of angles in [0, π). The wavelet at scale j and angle θ is defined by
ψj,θ(u) = 2 −2jψ(2−jrθu), u ∈ R2
Scales equal or larger than J are carried by a low-pass filter φJ .
The wavelet transform of an image x ∈ RN×N is a family of functions obtained by the convolution of x with discrete wavelets.4 Let Λ := {0, · · · , J − 1} × πL{0, · · · , L − 1} be an index set. The wavelet coefficients are
x ? ψj,θ(u) = ∑ v∈ΩN x(u− v)ψj,θ(v), u ∈ ΩN , (j, θ) ∈ Λ. (2)
The low-pass coefficients x ? φJ are defined similarly. 2The complex-valued representation Rx(γ, u) ∈ C is a function of (γ, u) in an index set Γ× ΩN . 3Such as Portilla & Simoncelli (2000); Gatys et al. (2015); Ustyuzhaninov et al. (2017); Zhang & Mallat (2021) 4The continuous wavelets are discretized with periodic boundary conditions on the spatial grid ΩN .

heading: 2.2.2 WAVELET PHASE HARMONICS AND THE PS MODEL
text: To model natural textures, it has been shown (Portilla & Simoncelli, 2000; Zhang & Mallat, 2021) that it is crucial to capture statistical dependencies between wavelet coefficients across scales. This can be achieved by using a wavelet phase harmonic representation, which is defined by the composition of a linear wavelet transform of x, and a non-linear phase harmonic transform.
In Mallat et al. (2020), the authors introduce the phase harmonics to adjust the phase of a complex number z ∈ C. More precisely, the phase harmonics {[z]k}k∈Z of a complex number z ∈ C are defined by multiplying its phase ϕ(z) of z by integers k, while keeping the modulus constant, i.e. ∀ k ∈ Z, [z]k := |z|eikϕ(z). The wavelet phase harmonic representation (WPH) is then defined by
RWPHx(γ, u) = [x ? ψj,θ(u)] k − µγ , γ = (j, θ, k) ∈ Γ = Λ× Z, (3)
where µγ is defined as the spatial average of RWPHx̄(γ, ·). It is shown in Zhang & Mallat (2021) that the PS model can be regarded as a low-order wavelet phase harmonics covariance model, which considers only a restricted number of pairs (k, k′) (see Appendix A for more details). In the next section, we shall use a dual representation of the phase harmonic operator to define a covariance model to capture high-order phase harmonics.

heading: 3 GENERALIZED RECTIFIER WAVELET COVARIANCE MODEL
text: In the previous section, we presented a class of models, built from the wavelet phase harmonic representation. A dual representation of the phase harmonic operator [·]k can be defined via a generalized rectified linear unit, that we review in Section 3.1. We then discuss in Section 3.2 how to define an appropriate index set of Γ for gray-scale textures. Section 3.3 extends the model to color textures.

heading: 3.1 FROM PHASE HARMONICS TO THE GENERALIZED RECTIFIER
text: The generalized rectified linear unit of a complex number z, with a phase shifted by α ∈ [0, 2π], is defined by
ρα(z) = ρ(Real(eiαz)), (4)
where ρ is a rectified linear unit, i.e. for any t ∈ R, ρ(t) := max(0, t). In Mallat et al. (2020), it is shown that applying a Fourier transform on ρα(z) along the variable α results in the phase harmonics of z (up to some normalization constant). This suggests an alternative model, defined by coefficients of the form
RALPHAx(γ, u) = ρα(x ? ψj,θ(u))− µγ , γ = (j, θ, α), (5)
for γ ∈ Γ = Λ× [0, 2π], and, as before, µγ is defined as the spatial average of RALPHAx̄(γ, u) over u ∈ ΩN .
Relation with high-order phase harmonics Based on the duality between the phase harmonics k ∈ Z and the phase shift variable α ∈ [0, 2π], we now present the relation between CALPHA and the high-order phase harmonics in CWPH. The following result was first proved in Mallat et al. (2020).
Proposition 1 There exists a complex-valued sequence {ck}k∈Z such that for all (j, θ, α) ∈ Γ, (j′, θ′, α′) ∈ Γ, and all τ ∈ ΩN ,
CALPHAx((j, θ, α), (j′, θ′, α′), τ) = ∑
(k,k′)∈Z2 ckc ∗ kC WPHx((j, θ, k), (j′, θ′, k′), τ)ei(kα−k ′α′).
The proof is given in Appendix B. We remark that the sequence {ck}k∈Z is uniquely determined by the rectifier non-linearity ρ, and they are non-zero if k is even (Mallat et al., 2020). This result shows that for a suitable choice of (α, α′). the covariance statistics CALPHAx can implicitly capture CWPHx with a wide range of k and k′.
Relation with second order statistics Using a simple decomposition of wavelet coefficients into their positive, negative, real and imaginary parts, we can further show that the covariance statistics CALPHAx capture the classical second order statistics of wavelet coefficients, also used in the PS model (with phase harmonic coefficients k = k′ = 1).
Proposition 2 Let I = {0, π2 , π, 3π 2 }. Their exists a finite complex-valued sequence {wα,α′}(α,α′)∈I2 such that for all (j, θ) ∈ Λ, and all τ ∈ ΩN ,∑ (α,α′)∈I2 wα,α′C ALPHAx((j, θ, α), (j′, θ′, α′), τ) = 1 ΩN ∑ u∈ΩN ( x ? ψj,θ(u) )( x ? ψj′,θ′(u− τ) )∗ .
(6)
The proof is given in Appendix C. This shows that using only four α uniformly chosen between [0, 2π] is sufficient to capture second order statistics. Because the wavelet transform is an invertible linear operator (on its range space), computing the r.h.s of eq. (6) for all (j, θ, τ), as well as the low-pass coefficients carried out by ΦJ , is equivalent to computing the correlation matrix of x.
Relation with the RF model Setting aside the subtraction by the spatial mean µγ , the RF model can be viewed as a particular case of models defined by eq. (5). Indeed, the statistics of the RF model take the form of eq. (1), with
RRFx(f, u) = ρ(x ? ψf (u)),
where {ψf} being a family of multi-scale random filters. By writing ρα(x ? ψj,θ(u − τ)) = ρ(x ? Real(ψτj,θe
iα)(u)), with ψτj,θ denoting the translation of ψj,θ by τ , we see that the models are similar, the difference being that our models use wavelet-based filters instead of random ones.

heading: 3.2 DEFINING AN APPROPRIATE Υ
text: The choice of the covariance set Υ is of central importance in the definition of the model. Intuitively, a too small set of indices will induce a model that could miss important structural information about the texture that we want to synthesize. Conversely, if Υ contains too many indices, the syntheses can have good visual quality, but the statistics of the model may have a large variance, leading to the memorization of some patterns of the observation. There is a trade-off between these two aspects: one must capture enough information to get syntheses of good visual quality, but not much, so as not to reproduce parts of the original image. To illustrate this point, we shall study the model ALPHA defined with three different sets Υ : A smaller model ALPHAS with a limited amount of elements in Υ, an intermediate model ALPHAI, and a larger model ALPHAL.
To precisely define these models, let us note J := {0, · · · , J−1}, Θ := πL{0, · · · , L−1}, andAA = 2π A {0, · · · , A−1}. Let us also define the set T := {0}∪{2
j(cos(θ), sin(θ))}0≤j<J, θ∈ πL{0,··· ,2L−1}, from which the spatial shift shall be selected. Table 1 summarizes the conditions that all parameters have to satisfy to be contained in these sets. Additionally, these models include large scale information through the covariance of a low-pass filter, i.e. the spatial average of x ?φJ(·)x ?φJ(· − τ), for τ ∈ T. To count the size of Υ without redundancies, Appendix A.3 provides an upper bound on the non-redundant statistics in our models. This upper bound is used to count the number of statistics in our models. To keep this number from being too large, instead of taking all shifts in a square box, such as in Portilla & Simoncelli (2000), we choose to select only shifts of dyadic moduli, and with the same orientations as the wavelets.
ALPHAS vs. ALPHAI The small model ALPHAS is inspired from the PS model, as it only takes into account of the interactions between nearby scales (i.e. |j′ − j| ≤ 1), and the spatial shift correlations are only considered for (j, θ) = (j′, θ′). There are two notable differences in the statistics included in the ALPHAS and ALPHAI models. The first one is the range of scales being correlated. It has been shown in Zhang & Mallat (2021) that constraining correlation between a wider range of scales induces a better model for non-Gaussian stationary processes, and a better estimation of cosmological parameters from observed data (Allys et al., 2020). The second difference, which has a significant impact on the number of statistics (it increases the model size by a factor ∼10), is the number of spatial shifts in the correlations. In the ALPHAI model, spatially shifted correlations
are computed for all pairs of coefficient (γ, γ′). For both stationary textures and non-stationary images in gray-scale, shape and contours of salient structures and objects are better reproduced with ALPHAI, as illustrated in Figure 2. More examples are given in Appendix D.
ALPHAI vs. ALPHAL As we observe in Figure 2, the ALPHAI model, containing 4 times less coefficients than the ALPHAL, suffers less from memorization effects, while still capturing most of the geometric information in the images. This small loss of information can be partially explained by the frequency transposition property of the phase harmonics operator (Mallat et al., 2020), for compactly supported wavelets in the frequency domain, as detailed in Appendix E. In order to avoid this memorization effect, we shall, in the rest of the paper, consider only the intermediate model.

heading: 3.3 MODELLING COLOR INTERACTIONS
text: In order to generate color textures, the covariance model ALPHAI defined in Section 3.2 could be directly applied to each R, G and B color channel independently. However, it would not take into account the color coherence in the structures of the observation.
To capture color interactions in the observation image, we shall impose the covariance between the coefficients of eq. (5) for all indices in Υ and all color channels. More precisely, let x = {xc}c=1,2,3 be a color image, with the parameter c representing the color channel. The ALPHAC color model is defined by correlations between coefficients of the form:
RALPHACx(γ, u) = ρα(x c ? ψj,θ)(u)− µγ , γ = (j, θ, α, c). (7)
The set of indices is defined as ΥALPHAC := {(γ, γ′, τ) : ((j, θ, α), (j′, θ′, α′), τ) ∈ ΥALPHAI , (c, c′) ∈ {1, 2, 3}2}.
Reduced ALPHAC The model ALPHAC has a large size as it computes correlations between all coefficients for all color channels. This size can be significantly reduced by computing spatially shifted coefficients only for the same color channels (to capture their geometries). This reduced model contains three times less coefficients (∼113k), without significant reduction of the visual quality of syntheses, as detailed in Appendix F.

heading: 4 NUMERICAL RESULTS
text: In this section, we compare our intermediate model to the state-of-art models (PS, RF and VGG) on both gray-scale and color textures. We first specify the experimental setup. We then present the synthesis results of various examples, and discuss their quality through visual inspection. A quantitative evaluation of the quality of the syntheses, based on the synthesis loss of the VGG model, and proposed in Ustyuzhaninov et al. (2017), is discussed in Appendix G.

heading: 4.1 EXPERIMENTAL SETUP
text: For our experiments, we choose gray-scale and color textures with various visual structures.5 In the gray-scale examples, we also include a stationary turbulent field (vorticity), which is simulated from Navier-Stokes equations in two dimensions (Schneider et al., 2006). These examples all contain complex geometric structures that are hard to model by the classical PS model.
The texture images presented in this work have a size of N = 256, giving a number of pixels of ∼65k. For all the ALPHA models, we use Morlet wavelets with a maximal scale J = 5 and number of orientations L = 4. This choice differs from the PS model, which uses Simoncelli steerable wavelets. In Appendix H, we discuss the impact of these wavelets on our model. To draw the samples, we follow gradient-based sampling algorithms, suitable in high-dimensions (Bruna & Mallat, 2018), to minimize the objective function ‖Cx−Cx̄‖2, starting from an initial sample from a normal distribution. Similarly to Gatys et al. (2015); Ustyuzhaninov et al. (2017), we use the LBFGS algorithm (Nocedal, 1980) for the optimization of the objective. As in the VGG model (Gatys et al., 2015), we further apply a histogram matching procedure as post-processing after optimization. The details of the PS, RF and VGG models, as well as detailed specifications of our models, are given in Appendix A. More synthesis examples can be found in Appendix I.

heading: 4.2 RESULTS
text: In Figure 3, we present examples of syntheses from the ALPHAI (or ALPHAC), PS, RF and VGG models, for both gray-scale and color textures, as well as for non-stationary images. We observe that our model ALPHAI produces texture syntheses of similar visual quality to the RF and VGG models. It also significantly outperforms the PS model in terms of the visual quality, without introducing visible memorization effects. As the model PS uses the statistics closer to ALPHAS compared to ALPHAI, the performance of PS is somehow expected.
5The source for the presented textures are given in Appendix A.
Note that for the tiles example (the fifth row) in Figure 3, the VGG model produces less convincing textures, because the long-range correlations present in the image (aligned tiles) are not reproduced. To remedy this issue, it has been proposed in Berger & Memisevic (2017) to add spatial shifts to the correlations of the network feature maps. These shift statistics are similar to the parameter τ in our model. We also observe that, in the case of the sixth row example (flowers), all models fail to reproduce complex structures at object-level. Possible improvements of such models is further discussed in Section 5.
For non-stationary images, we find that certain image patches can be more or less memorized by the RF, VGG and ALPHAI models, as illustrated in the seventh row example. Understanding such memorization effect of non-stationary images is a subtle topic, as we find that in some binary images (x̄(u) ∈ {0, 1}), only the PS model can reproduce the observation, even though it has a much smaller number of statistics (the last row example). We find that this is related to the spatial correlation statistics in PS (non-zero τ ). By removing this constraint, x̄ is no longer always reproduced.6 The non-stationary nature also appear in some logo near the boundary of some textures (e.g. bottom left in the observation of the first and fourth rows). Although this logo is reproduced by RF, VGG, ALPHAI and ALPHAC, it is a very local phenomenon, as we do not find visible copies of the textures when there is no logo, and it is likely related to the way one addresses the boundary effect (see more in Appendix A.5).

heading: 5 DISCUSSION
text: In this work, we presented a new wavelet-based model for natural textures. This model incorporates a wide range of statistics, by computing covariances between rectified wavelet coefficients, at different scales, phases and positions. We showed that this model is able to capture and reproduce complex geometric structures, present in natural textures or physical processes, producing syntheses of similar quality to state-of-the-art models that use CNN-based representations.
Being defined with a wavelet family instead of multi-scale random filters, the proposed model uses less statistics than the RF model. For the gray-scale textures, our model has about 15 times less statistics, as it focuses on capturing the geometric structures present in images. Although our color model has a slightly larger number of statistics than VGG, the reduced color model, presented in Section 3.3 is three times smaller than ALPHAC, while achieving similar visual quality. It shows the potential to further reduce the size of the color model. In the PS model, a PCA on the color channels is performed (Vacher & Briand, 2021). The same idea could also be applied to our model.
Furthermore, there are examples where all the models may all fail to produce some geometric structures at object-level, as illustrated in Figure 3 (sixth row). In this situation, we still need to find more informative statistics. One may for example consider to incorporate a second layer of wavelet transform as in the wavelet scattering transform (Leonarduzzi et al., 2019; Sifre & Mallat, 2013). Another line of research is to introduce other kinds of losses (such as to encourage image smoothness) in order to improve the VGG model (Sendik & Cohen-Or, 2017; Liu et al., 2016; Zhou et al., 2018). These losses are complementary, and thus could also be applied to our models.
Finding a minimal set of statistics to define a texture model remains important because a large number of statistics can result in a high variance of the estimators, and the associated model may suffer from memorization effects. This is a problem because the aim of the model is to approximate the underlying distribution of the observation, and therefore produce diverse textures. In this regard, the mere visual evaluation of the synthetic textures can fail to take this aspect of the model into account. A qualitative evaluation, coherent with visual perception, remains an interesting line of research to study the optimal trade-off between quality and diversity (Ustyuzhaninov et al., 2017).
Reproducibility Statement In order to reproduce the experiments presented in this work, the main parameters to define the model can be found in Section 4.1. More details are given in Appendix A. Additionally, the code shall be made publicly available.
6Set the parameter ∆ = 0 in the PS model. See Appendix A for more details. This simple example suggests that sometimes it is very important to choose the right statistics to capture specific geometric structures.
Ethics Statement The authors acknowledge that no potential conflicts of interest, discrimination, bias, fairness concerns or research integrity issues have been raised during the completion of this work.

heading: A MODEL AND ALGORITHMIC SPECIFICATION
text: We provide additional information needed to reproduce the numerical results of the models (for both gray-scale and color textures) considered in this paper. First, we detail the models of PS, RF, VGG. We then give algorithmic parameters to obtain the synthesis images. For natural textures, we also propose a strategy to synthesize non-periodic images in our models. An image is non-periodic if a periodic extension of the image to the domain outside ΩN create discontinuities at the contour of ΩN .
Sources of textures Our natural texture examples were obtained from the following three sources: CNS NYU7, Textures.com8, Describable Textures Dataset model9 and the Github page of Berger & Memisevic (2017) 10.
7http://www.cns.nyu.edu/˜lcv/texture/ 8https://textures.com/ 9https://www.robots.ox.ac.uk/˜vgg/data/dtd/index.html
10https://github.com/guillaumebrg/texture_generation
A.1 MODEL PARAMETERS
We specify the model parameters to synthesize both gray-scale and color textures. We also discuss how to extend the RF and the VGG model, originally designed for color textures, to model gray-scale textures.
• PS: For both gray-scale and color models, we set the number of scales J = 5, and the number of orientations L = 8 for the Simoncelli steerable wavelets. The spatial shift τ = (τ1, τ2) ∈ ΩN is chosen to be in the range of max (τ1, τ2) ≤ ∆ = 4. Note that this is different to the model parameters reported for the default PS model (which is J = 4, L = 4,∆ = 3) as we find that it results in a larger set of statistics and better visual quality. The synthesis results of this model can be reproduced by a Matlab software.11
• RF: For gray-scale textures, we consider J × L random convolutional filters (ψf )1≤f≤JL. Let f = (j, `), the index j representing the scale of each filter, whose size is (Wj ,Wj) for 1 ≤ j ≤ J . For Γ = {f = (j, `)|j ≤ J, ` ≤ L}, the representation is
RRFx(γ, u) = ρ(x ? ψj,`(u)), γ = (j, `) ∈ Γ.
For a color image x = {xc}1≤c≤3, we use 3 × J × L random convolutional filters. The representation of x is
RRFx(γ, u) = ρ ( 3∑ c=1 xc ? ψc,j,`(u) ) , γ = (j, `) ∈ Γ.
The correlations CRFx are defined for all pairs (γ, γ′) ∈ Γ× Γ. In both the gray and color cases, it results in a correlation matrix CRF with J2L2 statistics. Following the default setting of the RF model, we set J = 8 and L = 128 for filters whose sizes are W1 = 3,W2 = 5,W3 = 7,W4 = 11,W5 = 15,W6 = 23,W7 = 37,W8 = 55. Each filter ψj,` or ψc,j,` is generated randomly according to GlorotUniform in the software Lasagne. 12
• VGG. For a color image x = {xc}, the VGG model computes a correlation matrix CVGGx between the features maps within different layers of a pre-trained CNN network. To adapt this model to gray-scale textures, we shall add one input layer which converts a gray-scale image y into a color image, by setting xc = y for each color channel c. This allows one to use the same CVGGx to compute the gradient of the VGG loss with respect to y, and therefore to synthesise a gray-scale texture. For both gray-scale and color textures, we use only five layers ’conv1 1’, ’pool1’, ’pool2’, ’pool3’, ’pool4’, as proposed in the original work.
• ALPHA. See the main text.
A.2 ALGORITHMIC PARAMETERS
We specify the optimization parameters used to synthesize both gray-scale and color textures.
• PS: It utilizes iterative projections onto constraint sets to generate textures. We set the number of iterations to 200.
• RF: It uses the L-BFGS procedure13 with a memory size 20, and with a maximal number of iterations 2000. The initialization for each pixel value of a gray-image is Uniform between [−1, 1]. For the color image case, each RGB channel is initialized independently with a Uniform distribution between [−1, 1]. To address non-zero mean textures (i.e. E(X(u)) 6= 0), the empirical mean of x̄ is subtracted from the input x to compute the representation. It is added back to the output of the optimization to produce a synthesis.
11https://www.cns.nyu.edu/˜lcv/texture/ 12https://lasagne.readthedocs.io/ 13scipy.optimize.fmin_l_bfgs_b in Python
• VGG: It uses the L-BFGS procedure14 with a memory size of 20, and with a maximal number of iterations of 2000.For both gray-scale and color images, bounds constraints are used for the optimization. The initialization of each pixel value is the standard normal distribution (zero mean, unit variance). To address non-zero mean textures (i.e. E(X(u)) 6= 0), the VGG mean is subtracted from the input x of the representation15. It is added back to the output after the optimization to produce a synthesis (with an additional histogram matching post-processing).
• ALPHA: For all the models in ALPHA, we use the L-BFGS optimization algorithm with restarts. Starting from the standard normal distribution, with mean and standard deviation estimated from the observation, we use the L-BFGS procedure implemented in Pytorch. It runs for 500 iterations and then it is restarted with an initialization obtained from the previous L-BFGS result. This is repeated 10 times to obtain the synthesis (with an additional histogram matching post-processing).
A.3 NUMBER OF STATISTICS IN THE ALPHA MODELS
In this section, we detail the number of statistics in the different ALPHA models. We begin by giving the formula for each model (note that we do not include the low-pass statistics, which numbers are negligible). To (partially) avoid redundancy in the coefficients, for all models, we compute only the correlations for indices in Υ such that j2 ≥ j1. This gives us the following formulas for the number of statistics:
• #(ALPHAS) = (2J − 1)|Θ4|2|A4|+ J |Θ4||A4||T| • #(ALPHAI) = 12J(J + 1)|Θ4| 2|A4||T|
• #(ALPHAL) = 12J(J + 1)|Θ4| 2|A4|2|T|
• #(ALPHAC) = 92J(J + 1)|Θ4| 2|A4||T|
Note that, for all ALPHA models, we also compute first order statistics µγ , i.e. the spatial averages of RALPHAx̄(γ, u). There are J |Θ4||A4| statistics of this sort in every model, which is negligible with respect to the total number of second order statistics.
Note also that there are still some redundancies in these statistics, as for j = j′, all correlations for θ, θ′, α, α′, τ are counted twice. The number of such statistics is16, for all model, superior to the number of first order statistics, which shows that our formula is in fact an upper bound for the exact number of statistics.
A.4 NUMBER OF STATISTICS IN THE PS MODELS
The PS model can be interpreted as a particular case of the WPH covariance model as it contains the following three key categories of statistics:
• Raw coefficient correlations (k = k′ = 1): they capture 2nd order statistics of a stationary process X , i.e. the correlations between X(u) and X(u′) for (u, u′) ∈ Ω2N . • Coefficient magnitude statistics (k = k′ = 0): they capture information of X beyond the 2nd order statistics. Very often, nearby scales and angles are considered in the model, such as j′ ≈ j and θ′ ≈ θ. • Cross-scale phase statistics (k = 1, k′ = 2): they capture local phase alignments of the wavelet coefficients at nearby scales j′ = j + 1, which are complementary to the magnitude correlations.
The major issue to count the number of statistics of this model is to avoid double counting the statistics which are the same. This is mostly due to the symmetries of the covariance matrices. We
14scipy.optimize.minimize in Python 15For the color image whose pixel value is between zero and one, the mean of BGR is 0.40760392, 0.45795686, 0.48501961. For the gray-scale, we simply take the average of the BGR mean. 16The number of such moments is of the order of J |Θ|2|A4| for the small model, J |Θ|2|A4|2 for the intermediate and large models, and 9J |Θ|2|A4|2 for the color model.
next detail how obtain the number of statistics for both the gray-scale and color model, by following the work of Portilla & Simoncelli (2000) and Vacher & Briand (2021).
As before, we assume the number of wavelet scales is J , the number of wavelet orientations is L, and the spatial shift range is within a square of size (2∆+1)×(2∆+1). Let us denoteNa = 2∆+1. We next describe in detail the number of statistics of each category counted in our paper,
A.4.1 PS IN GRAY-SCALE
•Marginal statistics of x: 6. They include mean, variance, skewness, etc. •Marginal statistics of wavelet coefficients: 2(J + 1) + 1. • Auto-correlation of wavelet coefficients (raw coefficient correlations): (J + 1)(N2a + 1)/2. • Auto-correlation of magnitude of wavelet coefficients (coefficient magnitude statistics): JL(N2a + 1)/2 + JL(L− 1)/2 + (J − 1)L2. •Mean of magnitude of wavelet coefficients: JL+ 2. This is not counted in the paper of Portilla & Simoncelli (2000), but it used in the Matlab software.
• Cross-correlation of phase of wavelet coefficients (cross-scale phase statistics): 2(J−1)L2 +JL2. The extra JL2 coefficients are cross-correlation of real sub-band cousin coefficients, which are not counted in the paper, but used in the Matlab software.
To compare with the model size of the original work of Portilla & Simoncelli (2000), one can check that the sum of these number is 792 when J = 4, L = 4,∆ = 3 (Na = 7). If we do not count the coefficients which are only counted in the Matlab software (JL + 2 + JL2 = 82), it results in 710 as reported in the original paper.
A.4.2 PS IN COLOR
•Marginal statistics of x and PCA transform of x: 6× 3 + 3× 4 = 30. •Marginal statistics of wavelet coefficients: 6(J + 1) + 9. • Auto-correlation of wavelet coefficients (raw coefficient correlations): 3(J + 2)(N2a + 1)/2. This is called Central autoCorr of the PCA bands in the Matalb software (which includes lowband).
• Auto-correlation of magnitude of wavelet coefficients (coefficient magnitude statistics): 3JL(N2a + 1)/2 + J(3L)(3L− 1)/2 + (J − 1)(3L)2. •Mean of magnitude of wavelet coefficients: 3(JL+ 2). • Cross-correlation of phase of wavelet coefficients (cross-scale phase statistics): J(3L)2 + (J − 1)(3L)(6L).
A.5 NON PERIODIC BOUNDARIES IN NATURAL IMAGES
The convolution operation in the wavelet transform (equation 2) is performed using the Fast Fourier Transform. Additionally, recall from Section 2.1 that spatial shifts are defined with periodic boundary conditions. This implies periodicity of the input image x. However, natural texture images are not periodic, so one needs to adapt the computation of coefficients to take into account possible border effects. To that end, instead of averaging over all u ∈ ΩN as in eq. (1), each correlation coefficient is averaged over a sub-window inside ΩN , which size depends on the scales of the coefficients being correlated. More precisely, let γ = (j, θ, α) and γ′ = (j′, θ′, α′). Note jm := max(j, j′). We define Ωjm := {u = (u1, u2) ∈ ΩN : 2jm ≤ ui < N − 2jm , i = 1, 2}. Then, for non periodic images, we compute
CALPHAx(γ, γ′, τ) = 1 |Ωjm | ∑ u∈ΩN 1Ωjm (u)1Ωjm (u− τ)R ALPHAx(γ, u)RALPHAx(γ′, u− τ),
(8)
where the spatial shifts are defined periodically. Note that the spatial averages µγ and µγ′ are also performed on Ωjm .

heading: B PROOF OF PROPOSITION 1
text: Using the fact that
ρ(Real(zeiα)) = ρ(Real(|z|ei(ϕ(z)+α))) = |z|ρ(cos(α+ ϕ(z)),
and computing the Fourier coefficients of the 2π-periodic function ρα(z) in the variable α, we obtain
F(ρα(z))(k) := 1
2π ∫ [0,2π] ρα(z)e −ikαdα
= |z| 1 2π ∫ [0,2π] ρ(cos(α+ ϕ(z))e−ikαdα
= |z|eikϕ(z)ck = [z]kck,
where ck is the Fourier transform of h(.) := ρ(cos(.)) at the frequency k. The function α 7→ ρα(z) being periodic in α, we have its decomposition in Fourier series
ρα(z) = ∑ k∈Z F(ρα(z))(k)eikα
= ∑ k∈Z ck[z] keikα.
We can then write, for any z, z′ ∈ C, and α, α′ ∈ [0, 2π],
ρα(z)ρα′(z ′)∗ = ∑ k,k′∈Z2 ckc ∗ k[z] k[z′]−k ′ ei(kα−k ′α′).
Replacing z and z′ by any two wavelet coefficients x ?ψj,θ(u) and x ?ψj′,θ′(u− τ), we thus obtain the relation in Proposition 1.

heading: C PROOF OF PROPOSITION 2
text: Let z ∈ C, and recall from eq. (4), that ρα(z) = ρ(Real(zeiα)). Note that we have the following relation
z = ρ0(z)− ρπ(z)− i(ρπ2 (z)− ρ 3π2 (z)). (9)
We can then write zz′∗ = ( ρ0(z)− ρπ(z)− i(ρπ2 (z)− ρ 3π2 (z)) )( ρ0(z ′)− ρπ(z′)− i(ρπ2 (z ′)− ρ 3π 2 (z′)) ) =
∑ α,α′∈I2 wα,α′ρα(z)ρα′(z ′),
with I = {0, π2 , π, 3π 2 }. Replacing z with x ? ψj,θ(u), z ′ with x ? ψj′,θ′(u − τ), and injecting this relation in eq. (1) gives us the desired result.
D SUPPLEMENTARY RESULTS FOR THE GRAY-SCALE ALPHA MODELS
Here, we present further visual comparison between the different ALPHA models for gray-scale images, to illustrate the trade-off between quality and diversity.
E RELATION BETWEEN ALPHAL AND ALPHAI
Here, we informally explain why, under some conditions on the wavelet family, setting α′ ∈ {0} in the ALPHA models should not lose too much (but still some) information captured by the statistics.
First, remark that the simple linear relation( ∑ α∈A4 ρα(z)e ikα ) ρ0(z) = ∑ α∈A4 ρα(z)ρ0(z)e ikα (9)
tells us that computing all correlations for α ∈ A4 and α′ = 0 gives us at least the information contained in the r.h.s. of eq. (9).
Furthermore, if α ∈ A4 = {0, · · · , 3π4 }, we make the following approximation∑ α∈A4 ρα(z)e ikα ' ∫ [0,2π] ρα(z)e −ikαdα = 2πck[z] k.
Recall also from the proof of Proposition 2, thatF(ρα(z))(k) = [z]kck, where the Fourier transform is taken along the variable α. Therefore,
ρ0(z) = ∑ k∈Z ck[z] k.
Therefore, ( ∑ α∈A4 ρα(z)e ikα ) ρ0(z) ' ( ∫ [0,2π] ρα(z)e −ikαdα ) ρ0(z)
= 2πck[z] k ( ∑ k′∈Z ck′ [z ′]k ′)
= 2πck ∑ k′∈Z ck′ [z] k[z′]k ′ .
Then, replacing z and z′ with wavelet coefficients, we get
∑ α∈AA e−ikαCALPHAx((j, θ, α), (j′, θ′, 0), τ) ' 2πck ∑ k′∈Z ck′C WPHx((j, θ, k), (j′, θ′, k′), τ).
Using Plancherel’s theorem, we can write that
CWPHx((j, θ, k), (j′, θ′, k′), τ) = 1 |ΩN | ∑
ω∈ 2πN ΩN
F([x ? ψλ]k)(ω)F([x ? tτψλ]−k ′ )(ω),
where the Fourier transform of an image x is defined by F(x)(ω) := ∑ u∈ΩN x(u)e
−iωu, and tτ denotes the translation by τ , i.e. tτf(·) = f(· − τ). Now, suppose that the wavelets ψα have disjoint compact frequency support, in balls Bλ(2−jC ′), where λ = 2−jr−θξ0, and ξ0 is the central frequency of the mother wavelet ψ (cf. ?). Suppose also that frequency transposition property of the phase harmonics operator (cf. Mallat et al. (2020)) is such that [x ? ψλ]k has (approximately) frequency support in Bkλ(k2−jC ′). Then, for all λ, λ′, and all k ∈ Z, there exists only one k∗ such that the frequency supports of [x ? ψλ]k and [x ? ψλ′ ]−k ∗ are not disjoint, i.e. such that CWPHx((j, θ, k), (j′, θ′, k′), τ) 6= 0. This tells us that∑ α∈AA e−ikαCALPHAx((j, θ, α), (j′, θ′, 0), τ) ' ckck∗CWPHx((j, θ, k), (j′, θ′, k∗), τ).
Thus, computing all correlations for α ∈ A4, and α′ = 0 gives us (approximately) all the information contained in WPH coefficients for any pair k, k′.
This result lies on several approximations, and strong assumptions about the wavelets, which are not fully met in practice. For this reason, setting α′ = 0 instead of α′ ∈ A4 effectively reduced the amount of information captured by the statistics, and therefore increases the diversity of the model. However, as we observe in Section 3.2, there is not too much information lost, and the resulting model still captures most of the important geometric structures in texture images.

heading: F REDUCED COLOR MODEL
text: One can reduce the number of statistics in the color model by selection the spatial shift parameter τ to be non-zero only for correlations between the same color channels. More precisely, it is defined by the following index set: Υ := {(γ, γ′, τ) : ((j, θ, α), (j′, θ′, α′), τ) ∈ ΥALPHAI , c = c′ ∈ {1, 2, 3}} ∪ {(γ, γ′, 0) : ((j, θ, α), (j′, θ′, α′), 0) ∈ ΥALPHAI , (c, c′) ∈ {1, 2, 3}2}. This gives a model of size ∼ 113k, with little degradation of the visual quality, as shown in Figure 5.
G VGG SCORE
In Ustyuzhaninov et al. (2017), the authors proposed to use the synthesis loss of the VGG model to evaluate the quality of syntheses from any model. The goal is to define a quantitative, and more objective evaluation method than mere visual inspection. Since the VGG model produces syntheses almost indistinguishable from real textures, it is natural to consider its loss to asses the quality of a synthesis. We computed this loss for the first two examples of Figure 3 (radishes and cherries), and the frist two examples of Figure 9 (gravel and Turbulence flow). Note however that this loss is not exactly the same as the one used in Ustyuzhaninov et al. (2017), as the layers selected to compute
the loss are different. In this work, we chose to use the layers suggested in Gatys et al. (2015), (i.e. ’conv1 1’, ’pool1’, ’pool2’, ’pool3’, and ’pool4’ of the VGG-19 network (Simonyan & Zisserman, 2014)), and compute the relative VGG loss17.
We notice that this score is not always consistent with visual inspection, as there are texture examples and models for which the syntheses do not look much like the observation image, yet produce a small VGG loss (see e.g. the first and last rows of Figure 3, the RF model syntheses have the smallest loss). It should also be noted that the VGG loss reported on the VGG syntheses is not the synthesis loss after optimization, as a histogram matching (HM) procedure is performed as post-processing after optimization. We observed that the VGG loss of the syntheses from the VGG model after HM was considerably higher than the one for syntheses before it, while being visually very similar as illustrated in Figure 6. These observations suggest that the VGG score suffers from instabilities after reaching a certain level (that is, if the VGG loss is small enough, small perturbations of the values of the image pixels might have a strong impact on the loss).
H INFLUENCE OF THE CHOICE OF THE WAVELET TRANSFORM
H.1 INFLUENCE OF THE WAVELET FAMILY
In Section 3.2, we illustrated the importance of the set of indices Υ that define the wavelet coefficients being correlated. Another important role is played by the choice of the wavelets used in equation 2. As illustrated in Figure 7, this choice can have a visible impact on the quality of the textures. We observe that, while on the first example, the coherence of the structures appear similar for the three wavelet families, the second example shows that the wavelets used in Portilla & Simoncelli (2000) are less efficient in reproducing the contours of the objects (pebbles). While in our experiments, we chose to use the classical Morlet wavelets, an optimal choice for the wavelet family remains an open problem.
H.2 INFLUENCE OF SCALE PARAMETER
Recall from Section 2.2.1, the wavelet transform of an image x is defined by
{x ? ψj,θ, x ? φJ}0≤j<J, θ∈ πL{0,··· ,L−1}.
17Using the code from https://github.com/ivust/random-texture-synthesis/blob/ master/vgg_loss.py (function style loss relative).
The maximal scale parameter J also plays an important role in the definition of the wavelet transform. It determines the scales of the structures being captured by the transform. If this parameter is too small, large structures in the observation image might not be captured and reproduced in the model syntheses. Conversely, if J is too large, then the large scale statistics may have a high variance, inducing a memorization effect in the syntheses. Figure 8 illustrates this point on two examples from Section 4.2. By setting J = 4 (i.e. the maximal range of structures captured by the wavelets is of size 24 = 16), we observe on the first example that the larger structures (bubbles) are not well reproduced. When J is set to 6, the observation is almost identically reproduced by the synthesis. Similarly on the second examples, several parts of the synthesis appear very similar to ones in the observation. We found that a suitable trade-off consists in setting J = 5 for images of size N = 256.
I SUPPLEMENTARY RESULTS OF ALPHAI/ALPHAC VS. PS, RF, VGG
In Figure 9 we present additional syntheses on various examples, from the PS, ALPHAC, RF and VGG models. These examples can be viewed as random (Turbulence flow, tree bark, porous stone), structured (gravel, paisley pattern, tree leaf, school text), or inhomogeneous (crafted pattern of third row, but also the porous stone).
We see that ALPHAL again significantly improves the visual quality of the PS model on gray-scale textures such as the gravel and turbulence flow. The visual quality of RF and VGG seems also worse on Turbulence flow compared to ALPHAL. In some examples such as tree leaf, we find the synthesis of all the models are similar. On the inhomogeneous porous stone, non of the models give satisfying visual results.
Similarly, in Figure 10 are presented supplementary syntheses form the color models, for structured images (radishes, bubbles, flowers), quasi-periodic images (scales, honeycomb, bricks), and nonstationary images (feathers). As previously observed, for highly structured quasi-periodic images such as the bricks example, the VGG model fails to capture long-range correlation, which can be solved using the method of Berger & Memisevic (2017). Syntheses of non-stationary images exhibit memorization effects, as previously observed.

Review summaries:
This work proposed a texture synthesis framework using the rectified wavelet coefficients. The paper claims that the proposed method cand achieve similar quality with the VGG feature based method (Gatys et al. 22015) and random filter based method (RF, Ustyuzhaninov et al 2017) and gets better quality than PS (Portilla & Simoncelli 2000, while requiring less number of statistics than RF. 

This paper presents a new image representation model based on wavelets and non-linear rectifiers that allows to synthesize complex geometric textures with a better visual quality than previous wavelet-based models.  The main interest of the paper is the usage of the mathematical results from Mallat et al 2020 and Zhang and Mallat 2021 on wavelet phase harmonics to image texture synthesis.   They also show that the PS model (Portilla & Simoncelli 2000) is a particular case of the wavelet phase harmonics-based (WPH) model. Both methods underline the importance of statistical dependencies between wavelet coefficients across scales. The specific choice of the parameters available for WPH is important to balance between reproducing good structural information and memorizing patterns from the observation.