-M Morel





Chapitre  : Peut-on faire de la synthèse à base de grandes textures ?

). Contrairement aux méthodes de synthèse par l'exemple classique qui utilisent des images d'entrée très petites le but est de synthétiser de très grandes images de texture. On étudie ce problème avec  • L. Raad, A. Desolneux and J.-M. Morel. Midway blending as a new technique for patch quilting. To be submitted.

• L. Raad, A. Desolneux and J

Introduction

Motivation

La synthèse de texture par l'exemple a pour but de générer à partir d'un échantillon de texture de nouvelles images qui sont perceptuellement équivalentes à celle de départ. Les méthodes peuvent se regrouper en deux catégories : les méthodes paramétriques et les méthodes non-paramétriques à base de patchs. Le premier groupe a pour but de caractériser une image de texture à partir d'un ensemble de statistiques qui définissent un processus stochastique sous-jacent. Les résultats visuels de ces méthodes sont satisfaisants, mais seulement pour un groupe réduit de types de texture. La synthèse pour des images de textures ayant des structures très contrastées peut échouer. La deuxième catégorie d'algorithme découpe, puis recolle de manière consistente des voisinages locaux de l'image de départ pour générer de nouvelles configurations plausibles de ces voisinages (ou patchs). Les résultats visuels de ces méthodes sont impressionnants. Néanmoins, on observe souvent des répétitions verbatim de grandes parties de l'image d'entrée qui du coup peuvent être reproduites plusieurs fois. De plus, ces algorithmes peuvent diverger, reproduisant de fac ¸on itérative une partie de l'image de l'entrée en négligeant le reste.

Dans cette thèse on présente une approche combinant des idées des méthodes paramétriques et des méthodes non paramétriques à base de patchs. On l'appelle synthèse localement Gaussienne et elle est présentée dans le Chapitre . On préserve dans cette nouvelle méthode les aspects positifs de chaque approche : la capacité d'innover des méthodes paramétriques, et la capacité de générer des textures fortement structurées des méthodes non paramétriques à base de patchs. Pour ce faire, on construit un modèle Gaussien multidimensionnel des auto-similarités d'une image de texture. Une nouvelle image est générée patch par patch, o ù pour chaque patch une loi normale multidimensionnelle est estimée à partir des patchs de l'image d'entrée qui lui ressemble le plus, et ensuite échantillonné à partir de cette loi. Les images synthétisées sont donc différentes point à point de la texture de départ. Ainsi, on obtient des résultats qui sont visuellement supérieurs à ceux obtenus avec les méthodes paramétriques et qui sont comparables à ceux obtenus avec les méthodes nonparamétriques à base de patchs tout en utilisant une paramétrisation locale de l'image. La thèse s'attache aussi à résoudre une autre difficulté des méthodes à base de patchs : le choix de la taille du patch. Afin de réduire significativement cette dépendance, on propose une extension multi-échelle de la méthode au Chapitre .

Les méthodes à bases de patchs supposent une étape de recollement. En effet, les patchs de l'image synthétisée se superposent entre eux, il faut donc gérer le recollement dans ces zones. La première approche qu'on a considérée consiste à prendre en compte cette contrainte de superposition dans la modélisation des patchs. Cette approche est décrite au Chapitre . Les expériences montrent que cela est satisfaisant pour des images de textures périodiques ou pseudo-périodiques  et qu'en conséquence l'étape de recollement peut être supprimée pour ces textures. Cependant, pour des images de textures plus complexes ce n'est pas le cas, ce qui nous a menée à proposer une nouvelle méthode de recollement inspirée du transport optimal qui est détaillée au Chapitre .

Cette thèse se termine avec une étude complète de l'état de l'art en génération d'images de textures naturelles dans le Chapitre , qui tient lieu aussi de conclusion. L'étude que nous présentons montre que, malgré les progrès considérables des méthodes de synthèse à base d'exemples proposées dans la vaste littérature, et même en les combinant astucieusement, celles-ci sont encore incapables d'émuler des textures complexes et non stationnaires.

Dans cette introduction on donne une vue d'ensemble de tous les chapitres de la thèse avec leur résultats principaux et conclusions. On termine avec la liste de publications en lien avec ces travaux.

Chapitre  : Modèles localement Gaussien

Les principales approches de la synthèse d'images de textures sont les méthodes paramétriques et les méthodes non-paramétriques à base de patch. Le premier modèle de texture caractérise l'image par un ensemble de statistiques. L'algorithme de synthèse associé estime ces statistiques à partir de l'échantillon de texture et génère une nouvelle image. La texture générée est en général moins précise que celle d'entrée. Le second modèle se résume à une démarche intelligente de "copiercoller" qui met ensemble des morceaux de la texture d'entrée. La texture générée peut présenter des répétitions verbatim de grandes parties de l'échantillon de texture. Au Chapitre  nous proposons de mélanger les deux approches mentionnées en utilisant un modèle de texture localement Gaussien dans l'espace des patchs. Cela permet de synthétiser des textures qui sont partout différentes de celle d'entrée tout en ayant une meilleur qualité visuelle par rapport aux résultats de méthodes paramétriques.

Inspiré de la méthode d'Efros et Freeman [] nous synthétisons une nouvelle texture patch par patch. A différence de [] on va supposer que chaque nouveau patch synthétisé suit une loi normale multidimensionnelle qui est apprise à partir d'un ensemble de patchs pris dans l'échantillon de texture. Cet ensemble est constitué des patchs les plus semblables, pour la norme L 2 , au patch qui est en cours de synthèse. Les modèles localement Gaussien ont été utilisés pour le débruitage d'image [] avec des résultats très satisfaisants. Pour la synthèse de texture, cette approche nous permets de maintenir la cohérence entre patchs par rapport à l'image d'entrée tout en créant de nouveaux patchs qui n'existent a priori pas dans l'échantillon de texture mais qui leur sont visuellement semblables.

Le modèle de patch Les lois normales multidimensionnelles impliquées sont définies par leur vecteur moyen µ (x,y) et leur matrice de covariance Σ (x,y) . Pour un patch donné p (x,y) u , sous sa forme vectorielle, et qui est de taille 1 × n 2 , les paramètres de la Gaussienne sont estimés à partir d'un ensemble de patchs que l'on dénote U = {p (x i ,y i ) u , i = 1, . . . , m}, o ù les p (x i ,y i ) u , i = 1, . . . , m sont les m patchs les plus semblables à p u (x, y) selon la norme L 2 . Le vecteur p(x,y) u défini comme dans l'équation (.) suit une loi normale N (µ (x,y) , Σ (x,y) ). On observe que les variances de ces modèles confirment un degré d'innovation raisonnable pour les vecteurs échantillonnés.

p(x,y)

u = 1 √ m -1 m i=1
a i (p (x i ,y i ) u µ (x,y) ) + µ (x,y) , a i ∼ N (0, 1), i = 1, . . . , m

 L'algorithme de synthèse localement Gaussienne On enchaîne sur la description d'un algorithme de synthèse qui génère une nouvelle image de texture en collant entre eux des patchs échantillonnés selon des lois normales multidimensionnelles (.). On note cet algorithme par LG.

Cette méthode est itérative : les patchs sont synthétisés en raster-scan (de haut en bas, de gauche à droite). Le but de chaque itération est de générer un nouveau patch p (x,y) w qui est partiellement défini sur une région qu'on appelle "région d'intersection" (voir Figure ). La partie connue du patch est utiliser pour définir l'ensemble U qui déterminera la loi normale du nouveau patch. Le patch p (x,y) w est généré comme défini dans (.). La dernière étape consiste à recoller le nouveau patch dans l'image de sortie. Pour cela la méthode utilisée est celle proposée dans []. En conclusion, l'algorithme de synthèse proposé génère une nouvelle image de texture w qui est visuellement semblable à la texture d'entrée et pourtant elle n'est pas composée des patchs existant dans l'image d'entrée. Par conséquent, la méthode proposée est capable de diminuer quelques uns des aspects négatifs des méthodes de synthèse dites paramétriques et non-paramétriques à base de patchs. L'algorithme reste cependant dépendant du choix de la taille du patch n et de m la taille de l'ensemble U , comme on peut observer dans la Figure  et qui devrait être ajuster selon l'échantillon detexture. Notre algorithme a une faible complexité comparé aux méthodes de débruitage à base de patch comme [], []. Une possibilité pour réduire la dépendance de la méthode à la taille du patch est de travailler avec une approche multi-échelle. Cette extension est présentée au Chapitre . Contrairement aux méthodes de synthèse de texture paramétriques et tout comme les méthodes non-paramétriques à base de patchs, notre algorithme n'est pas forcé à respecter les statistiques globales de l'image d'entrée. Cet aspect est également considéré au Chapitre . Un dernier point tout aussi important est traité aux Chapitres  et  et qui concerne le recollement d'un nouveau patch dans une partie de l'image déjà synthétisée. patchs. Notre but est que le modèle considéré prenne en compte l'étape de recollement des patchs. Aussi, nous proposons de généraliser l'algorithme de synthèse localement Gaussienne (LG en utilisant une approche multi-échelle o ù la texture est considérée comme étant localement Gaussienne aux différentes échelles. Cette approche est évaluée en utilisant une grande variété d'images de texture (réelles et synthétiques) et les résultats montrent l'efficacité de la métode proposée pour une grande gamme de texture.

Modèle conditionnel Gaussien

La première question qu'on s'est posé était si l'étape de recollement des patchs pouvait être supprimée et comment ? Dans le Chapitre , on a adopté l'approche d'Efros et Freeman [] o ù pour chaque nouveau patch synthétisé son bord, sur la région d'intersection, est redéfini. Dans ce chapitre on considère deux modèles qui sont contraints aux valeurs de la région d'intersection : le modèle Gaussien conditionnel (CLG) et le modèle Gaussien conditionnel régularisé (RCLG). Ces deux solutions montrent des résultats satisfaisants pour des images de texture périodiques et pseudo-périodiques mais elles ont tendance a échoué pour des textures plus complexes. On peut voir dans la

Généralisation multi-échelle

La deuxième question qu'on s'est posé par rapport aux méthodes à base de patchs était comment gérer leur forte dépendance à la taille du patch utilisé en particulier quand on travaille avec des images fortement structurées. En effet, les textures structurées ont en général des détails à différentes échelles qui ne peuvent pas être capturés avec une seule taille de patch. L'approche multi-échelle que nous envisageons peut être résumée en quelques lignes et qu'on dénote M SLG. La méthode commence à l'échelle la plus grossière o ù on utilise l'algorithme présenté au Chapitre  avec une différence : l'étape de recollement des patchs est remplacée par une combinaison linéaire des valeurs sur la région d'intersection. Aux échelles suivantes la synthèse est faite telle qu'à chaque échelle les détails fins sont rajoutés à l'image synthétisée à l'échelle précédente. On suppose que les détails suivent une loi normale multidimensionnelle. En résumé, une image de plus faible résolution est synthétisée et utilisée comme guide à laquelle on rajoute les détails des échelles restantes. Une description complète est fournie dans le Chapitre . Dans la Figure  on montre quelques résultats de synthèse multi-échelle.

Chapitre  : Recollement midway de patchs

Dans le Chapitre  on s'intéresse à la phase de recollement de patchs des méthodes de synthèse à base de patchs. Un des défis de cette étape est de recoller des patchs qui diffèrent considérablement sur leur bord commun. Les techniques utilisées pour la synthèse de texture ont tendance à échouer dans des cas de fort contraste et par conséquent on observe la transition entre les patchs (voir Figure ).

On propose une nouvelle méthode qui se base sur la théorie du transport optimal afin de fusionner deux patchs tout en ayant une transition lisse entre eux. Cette approche peut se décomposer en deux étapes. La première consiste à définir un bord commun aux deux patchs qui se trouve à michemin entre leur bords initiaux. Pour obtenir ce bord à mi-chemin deux déformations différentes sont appliquées à chaque patch individuellement. La deuxième étape consiste à propager graduellement ces déformations sur chacun des patchs afin de générer un changement graduel et lisse sur chacun des patchs.

Le signal à mi-chemin La première étape consiste à trouver un signal de dimension un qui se trouve à mi-chemin entre deux signaux de dimension un qu'on dénote s 1 et s 2 . Ces signaux, s 1 et s 2 , correspondent aux bords des patchs p 1 et p 2 sur lesquels on veut les recoller (voir Figure ). On considère s 1 et s 2 comme étant deux histogrammes qu'on veut transporter dans l'esprit de la méthode d'égalisation d'histogrammes midway []. Le signal s à mi-chemin est défini tel que son histogramme cumulé soit la moyenne harmonique des histogrammes cumulés de s 1 et s 2 . Pour passer des signaux s 1 et s 2 au signal s on applique deux déformations différentes φ 1 et φ 2 à s 1 et s 2 qu'on appellera les déformations midway.

Propagation des déformations midway Une fois que les signaux s 1 et s 2 coïncident il faut propager les déformations φ 1 et φ 2 afin que dans les patchs p 1 et p 2 il n'y ait pas de discontinuités entre s et le reste du patch respectivement. Pour cela on considère deux alternatives. La  première consiste à propager linéairement et graduellement les déformations midway φ 1 et φ 2 dans chaque patch. La deuxième consiste à définir des signaux intermédiaires (pour chaque patch) qui résultent d'une combinaison convexe entre l'inverse de l'histogramme cumulé de s et l'inverse de l'histogramme cumulé du signal provenant du bord qui délimite la zone de propagation dans chacun des patchs. La première méthode de propagation rend des résultats o ù les déformations appliquées sont très lisses. Cependant, des problèmes de normalisation sont évidents et résultent du fait de considérer les signaux comme des histogrammes (qu'il faut donc normaliser pour qu'ils aient la même masse). La deuxième méthode est exempte de ces problèmes de normalisation. Néanmoins, les déformations sont bien plus brusques et peuvent être quelque fois plus visibles que dans le cas linéaire.

Nous avons évalué en détail les méthodes proposées. Dans une première partie nous avons analysé le recollement midway en sélectionnant deux patchs de manière aléatoire dans une même image de texture. Nous avons comparé le recollement midway à deux autres méthodes : Poisson editing [] et le quilting d'Efros et Freeman []. Quelques exemples difficiles sont montrés dans la Figure . Dans une deuxième partie nous avons testé le recollement midway dans l'algorithme de synthèse d'Efros et Freeman [] et comparé les résultats à ceux de la méthode originale. Ces résultats nous ont permis de conclure que la méthode proposée est une possible alternative pour le recollement des patchs dans les méthodes de synthèse à base de patchs. les outils qui sont à notre disposition : les méthodes qui ressortent de la vaste littérature de la synthèse par l'exemple ainsi que nos contributions des Chapitres  et , et des combinaisons de ces méthodes.

Pour anticiper nos conclusions aucune des méthodes existantes n'est capable de synthétiser correctement de grandes images de texture. Néanmoins, nous explorons des combinaisons astucieuses de ces méthodes qui présentent des résultats intéressants mais pas parfaits. Les méthodes existantes montrent des résultats très satisfaisants quand elles sont appliquées à de tout petit bout des grandes texture d'entrée (voir Figure ). Nous observons que ces grandes textures sont très différentes des textures académiques utilisées jusqu'à présent. Au fait elles montrent des variations locales très prononcées, qui peuvent être vues comme différentes textures, même si leur aspect global est celui d'une seule texture (Figure ). Cette étude est une bonne expérience pour mettre en épreuve toutes les méthodes de synthèse à base d'exemple classiques.

Cette étude commence par tester différentes méthodes de synthèse de texture par l'exemple notamment celles d'Efros et Freeman [], Portilla et Simoncello [], Galerne et al. [] et l'approche multi-échelle présenté au Chapitre  appliquées à différentes images de texture naturelles et de très grande taille comme du bois, de la pierre et du métal. On a observé que les résultats de synthèse sur les images de bois n'étaient pas satisfaisants. Néanmoins quand on a considéré des sous-images prises dans ces grandes textures de bois et o ù aucune structure saillante apparaît, les résultats de synthèse sont bons. On peut rapidement conclure que la perception qu'on a de ces images comme étant une seule texture n'est pas "vraie". Un exemple est montré avec les Figures  et  FIGURE 9 -On montre quatre sous-images prise de l'image de la Figure 1.8 à la r ésolution originale. Ces sous-images sont de taille 500 × 500 pixels. On observe que chacune de ces quatre images repr ésentent des textures bien diff érentes appartenants à une seule "grande" texture.

jeu d'expériences consiste à combiner ces méthodes de synthèse, en particulier :

• La méthode d'Efros et Freeman [] combinée à la synthèse localement Gaussienne (LG + EF )

• La méthode de Portilla et Simoncelli [] combinée à la synthèse localement Gaussienne (LG + P S)

Les combinaisons précédentes améliorent significativement les résultats obtenus dans la première étape en particulier la combinaison LG+PS. Pourtant cette alternative n'est pas suffisante dans le cas des images de bois. Par la suite on montre différentes tentatives pour synthétiser ces images comme par exemple l'application de déformation géométriques en espérant varier l'aspect des structures fortement saillantes. Une astuce qui est souvent utile est de permuter de manière aléatoire des sous-parties de l'image d'entrée afin de décomposer ses détails fortement reconnaissables comme dans le cas du bois et de les disposer de fac ¸on aléatoire dans la nouvelle image. Une dernière expérience que nous avons considérée est de varier les basses fréquences de l'image d'entrée en utilisant la méthode de Galerne et al. [] tout en laissant inchangé les hautes fréquences de l'image d'entrée. Dans ces trois essais de synthèse nous observons des comportements intéressants pour les différentes textures. Pour illustration trois exemples sont montrés dans la

Liste de Publications

• L. Raad, A. Desolneux and J.-M. Morel. Locally Gaussian Exemplar Based Texture Synthesis.

In IEEE International Conference on Image Processing (ICIP), pages  -, . IEEE.

• L. Raad, A. Desolneux and J.-M. Morel. Conditional Gaussian Models for Texture Synthesis.

In International Conference on Scale Space and Variational Methods in Computer Vision (SSVM), pages  -, . Springer International Publishing.

• L. Raad, A. Desolneux and J.-M. Morel. Multiscale Exemplar Based Texture Synthesis by Locally Gaussian Models. In IberoAmerican Congress on Pattern Recognition (CIARP),pages  -, . Springer Iternational Publishing.

• L. Raad, A. Desolneux and J.-M. Morel. Conditional and multiscale locally Gaussian models for texture synthesis. In Journal of Mathematical Imaging and Vision (JMIV), pages  -, .

• L. Raad and B. Galerne. Efros and Freeman Image Quilting Algorithm for Texture Synthesis. Submitted to Image Processing On Line (IPOL).

 1 Introduction

Motivation

Exemplar based texture synthesis is the process of generating, from an input texture sample, new texture images that are perceptually equivalent to the input. There are roughly two main categories of algorithms: the statistics-based methods and the non-parametric patch-based methods. The statistics-based methods aim at characterizing a given texture sample by estimating a set of statistics which will define an underlying stochastic process. The new images will then be samples of this stochastic process, i.e. they will have the same statistics as the input sample. The question here is what would be the appropriate set of statistics to yield a correct synthesis for the wide variety of texture images? The results of statistical methods are satisfying but only on a small group of textures, and often fail when important structures are visible in the input. The non-parametric patch-based methods reorganize local neighborhoods from the input sample in a consistent way to create new texture images. These methods return impressive visual results. Nevertheless, they often yield verbatim copies of large parts of the input sample. Furthermore, they can diverge, starting to reproduce iteratively one part of the input sample and neglecting the rest of it, thus growing what experts call "garbage".

In this thesis we propose a technique combining ideas from the statistical based methods and from the non-parametric patch-based methods. We call it the locally Gaussian method and it is presented in Chapter . The method retains the positive aspects of both categories: the innovation capacity of the parametric methods and the ability to synthesize highly structured textures of the non-parametric methods. To this aim, the self-similarities of a given input texture are modeled with conditional multivariate Gaussian distributions in the patch space. A new image is generated patch-wise, where for each given patch a multivariate Gaussian model is inferred from its nearest neighbors in the patch space of the input sample, and hereafter sampled from this model. The synthesized textures are therefore everywhere different from the original. In general, the results obtained are visually superior to those obtained with statistical based methods, which can be explained by the use of a local parametric model instead of a global one. On the other hand, our results are comparable to the visual results obtained with the non-parametric patch-based methods. This dissertation addresses in Chapter  another weakness of patch-based methods. They are strongly dependent on the patch size, which has to be decided manually. It is therefore crucial to fix a correct patch size for each synthesis. Since texture images have, in general, details at different scales, we extend the method to a multiscale approach which reduces the strong dependency of the method on the patch size.

Patch based methods involve a stitching step. Indeed, the patches used for the synthesis process overlap each other. This overlap must be taken into account to avoid any transition artifact  from patch to patch. Our first attempt to deal with it was to consider directly the overlap constraints in the local parametric model. This is described in Chapter . The experiments show that for periodic and pseudo-periodic textures, considering these constraints in the parametrization is enough to avoid the stitching step. Nevertheless, for more complex textures it is not sufficient. This leads us to suggest a new stitching technique inspired by optimal transport and midway histogram equalization which is detailed in Chapter .

This thesis ends with an extensive analysis of the generation of several natural textures in the conclusive Chapter . This study shows that, in spite of remarkable progress for local textures, the methods proposed in the extensive literature of exemplar based texture synthesis are still limited when dealing with complex and non-stationary textures.

In this introduction, we give a brief overview of all the chapters and their main results and conclusions. We end it with a list of publications linked to this PhD.

Chapter : Locally Gaussian models

The main approaches to texture modeling are the statistics-based methods and the non-parametric patch-based methods. In the first model the texture is characterized by a sophisticated statistical signature. The associated sampling algorithm estimates this signature from the example and produces a genuinely different texture. This texture nevertheless often loses accuracy. The second model boils down to a clever "copy-paste" procedure, which stitches together verbatim copies of large regions of the example. In Chapter  we propose to involve a locally Gaussian texture model in the patch space. It permits to synthesize textures that are everywhere different from the original but with better quality than the purely statistical methods.

We model each texture patch by a multivariate Gaussian distribution learned from its similar patches. Inspired by [], we maintain the idea of searching for patches to stitch together in the original sample. However, instead of using the exact patch taken in the input texture, we sample the stitched patch from its Gaussian model. The locally Gaussian patch models have been proved very useful in image denoising []. Our approach permits to maintain the coherence between patches with respect to the input sample, while creating new patches that do not exist in the sample texture but are still perceptually equivalent to it.

The patch model

The multivariate Gaussian models involved are defined by their mean vector µ (x,y) and their covariance matrix Σ (x,y) . For a given patch p (x,y) u , of size n × n pixels, these parameters are estimated from the set of the m nearest patches U = {p y) , Σ (x,y) ). We observed that indeed these models have reasonable variances confirming that effectively the patches simulated have an acceptable degree of innovation.

(x i ,y i ) u , i = 1, . . . , m}. The sampled vector p(x,y) u defined in (.) follows the distribution N (µ (x,

p(x,y)

u = 1 √ m -1 m i=1 a i (p (x i ,y i ) u -µ (x,y) ) + µ (x,y) , a i ∼ N (0, 1), i = 1, . . . , m (.)
The locally Gaussian synthesis algorithm We follow by describing a synthesis algorithm that generates a new texture image by stitching together patches sampled from multivariate Gaussian distributions (.) in the input sample patch space. This method is iterative: the patches are synthesized in a raster-scan order (top to bottom and left to right). The goal of each iteration is to that is partially defined on a region called the overlap area (see Figure .). The known part of the patch will define the set of patches U from which its Gaussian model is inferred. The generated patch p (x,y) ũ is then sampled as defined in (.). The last step consists in stitching the patch into the output texture using the quilting method of []. This synthesis algorithm generates a texture that is perceptually equivalent to the sample texture yet not composed of patches existing in the input texture. Thus, this method reduces some of the drawbacks of the statistics-based and the patch-based methods. The algorithm remains dependent on the choice of the patch size n and of the number of nearest neighbours m as can be observed in Figure ., that may have to be adjusted for each texture sample. Our algorithm has low complexity, compared for instance with classic patch-based denoising algorithms [], []. An alternative to reduce the dependency of the method to the patch size is to work in a multiscale approach. This extension is presented in Chapter . Unlike the statistics-based algorithms, but like the other patch-based methods, our algorithm is not forced to respect the global statistics of the texture sample. This issue is also considered in Chapter . A last and not minor aspect that will be addressed in Chapter  and Chapter  is the way the patches are actually stitched together.

Vertical overlap

Chapter  is structured as follows. In Section ., we give an estimation method of the Gaussian model for each patch. An analysis of the model's variance is provided. In Section ., the locally Gaussian synthesis algorithm is described. Section . compares the results with state of the art texture-synthesis methods. The concluding section . also states the limitations of the method.

Chapter : Conditional Gaussian model: a multiscale algorithm

In this chapter, we extend the model of Chapter  to model texture self-similarity with conditional Gaussian distributions in the patch space. Our goal is to extend the use of stitching techniques. In continuation the algorithm is also generalized by a multiscale procedure, where texture patches are modeled at each scale as spatially variable Gaussian vectors. This approach is tested over several real and synthetic texture images and its results show the effectiveness of the proposed technique for a wide range of textures.



Conditional Gaussian models

The first question that was brought up was how to stitch together the patches. In Chapter , the Efros and Freeman's approach [] was adopted, where every new patch overlaps the previously synthesized one, and the overlapped parts are blended. In this chapter we consider two models constrained to respect the values of its overlapping zone: the Conditional Gaussian model (CLG) and the Regularized Conditional Gaussian model (RCLG). This solution shows satisfying results for periodic or pseudo-periodic textures but sometimes fails for more complex ones. Two texture examples are shown in . where the three models LG, CLG and RCLG are compared.

Multiscale generalization

The second and central question brought up by non-parametric patch-based methods is: how to handle the strong dependency of the method on the patch size, in particular when dealing with macro-textures. Indeed, macro-textures show information at different scales that cannot be captured with a unique patch size. The multiscale approach can be summarized in a few sentences. The method begins by a synthesis at the coarsest scale using the local Gaussian method presented in Chapter  where the stitching step is replaced by a simple av-



input

LG model CLG model RCLG model Chapter  is structured as follows. In Section . the local Gaussian (LG) model is described as well as the description of its synthesis method. In Section . two new patch models are introduced: the conditional local Gaussian model (CLG) and the regularized conditional local Gaussian model (RCLG). Section . presents the multiscale approach (MSLG). Section . shows several experiments: a comparison of the three patch models proposed, a comparison to the results of state of the art texture synthesis methods and the influence of the parameters involved in the multiscale texture synthesis algorithm. Conclusions are presented in Section ..

Chapter : Midway patch blending

The problem of stitching patches together is the subject of Chapter . One of the major challenges is stitching patches that differ considerably along the edges. Stitching methods used in texture synthesis usually fail in these cases, resulting in evident contrast changes between the patches (see Figure .).

We suggest a novel method based on optimal transport theory that manages to merge both patches producing a smooth transition between them. This approach works in two steps: first by For both methods we used the same seed patch for each texture sample. defining a "midway" common restriction for both patches on their common stitchig edge. Midway also defines a different deformation along the edge for each patch. Then this deforamtion is propagated into both patches to produce a smooth transition.

The midway signal

The first step can be described as finding a common one dimensional signal which is midway, up to a minimal deformation, between two one dimensional signals denoted s 1 and s 2 (see Figure .). The signals under construction are nothing but the restriction of the patches p 1 and p 2 to the blending edge. The two signals s 1 and s 2 are handled as two histograms on which we apply a geometric deformation derived from the midway histogram equalization [].

The common intermediary signal s is thus defined as the underlying histogram of the harmonic mean of the cumulative functions of s 1 and s 2 .

Propagating the midway signal For the propagation step two alternatives are considered. The first one consists in linearly propagating the deformations. The second one considers convex combinations of two inverse cumulative signals: the midway signal and the signal on the curve up to which the deformation is propagated. The first alternative yields smoother results but is problematic regarding the normalization step of the signals (to impose the same mass). The second option avoids the normalization issues yet the propagation is more brutal and sometimes more  conspicuous than the linear alternative. We thoroughly evaluated the proposed method first by trying it directly on randomly selected patches and comparing it to two stitching techniques: the Poisson editing method [] and the Efros and Freeman method []. A few challenging examples examples are shown in Figure .. We then tried this method as a new stitching approach in the Efros and Freeman synthesis methods []. Based on the results, the proposed method proved to be a possible alternative to perform patch stitching in non-parametric patch-based methods. Chapter  is organized as follows. Section . recalls the midway equalization algorithm []. In section . the midway blending technique is introduced. At first a description on a simple blending curve is given. It is followed by a generalization of the method. A description of the algorithm is provided as well as a discussion on the extension to a color version. As a final contribution of this section we describe how to integrate this method into a simple and efficient texture synthesis method. In section . experimental results are organized in two parts. The first one compares the midway blending method to [] and []. The second part compares synthesis results of [] to the same synthesis method using the midway blending. Conclusions are given in section ..

Chapter : Can we emulate large textures?

In this conclusive chapter we aim at analyzing the limits of exemplar-based texture methods when using large natural texture images (see Figure .), a new problem which we might call Big Exemplar Based texture synthesis. Unlike exemplar-based methods that use small texture samples the goal is to synthesize large textures. We shall explore this problem with the tools at reach, namely the best methods from the extensive literature including our contributions of the previous chapters, and will not hesitate to combine them. To anticipate our conclusions, none of the existing methods manage by itself to satisfactorily emulate "big textures". Yet, we'll find that a clever combination of these methods yields interesting-if not perfect-results. Interestingly, the existing methods provide very good results on small crops (see Figure .) of "big textures" images. But we discovered with We began by trying some of the exemplar based texture synthesis approaches (Efros and Freeman [], Portilla and Simoncelli [], Galerne et al. [] and the multiscale locally Gaussian method presented in Chapter ) on different types of natural large texture images such as: wood, stone and metal. We observed that the synthesis results for the wood images were not satisfying. However for local crops of these textures where no strong global organization were evident the results obtained are impressive. One concludes obviously that, in spite of their perceptual unity, large textured images are more than a texture in the usual academic sense, as can be seen for the wood texture example shown in • The Efros and Freeman method [] combined to the locally Gaussian method (LG+EF)



• The Portilla and Simoncelli method [] combined to the locally Gaussian method (LG+PS) Combining the methods yields better visual results in particular for the combination PS + LG. Still this combination will not make it for wood textures, for example. We then show other emulation attempts, like for example the application of geometric deformations to the inputs aiming at changing the aspect of the salient structures. An often successful trick is randomly switching sub-parts of the input to partition the salient structures and reorganize them in a way that visually mitigates the repeated patterns. As a last experiment we tried to vary the low frequency information using the method proposed by Galerne et al. in [] while leaving unchanged the high frequencies with the same purpose. In all three emulation attempts we observed interesting behaviours for some of the textures. In Chapter  is organized as follows. Section . describes the large natural texture images we worked on. In section . a brief analysis of some of the state of the art exemplar-based texture synthesis methods is provided listing the pros and cons of each of them. In section . a first set of experiments are showed on small crops of the big textures. All of the exemplar-based methods of section . are tested on these images. They are tested separately and then some of them are combined. The results on these images are very satisfying. The same experiments are performed on the big texture themselves to try to emulate their structures at low resolution. The results are not very satisfying since the typical drawbacks of exemplar based methods arise (verbatim copies and garbage growing). This brings us to test other alternatives presented in Sections ., . and .. We try the application of one and two dimensional anamorphoses to the inputs, a random combination of slices of the input texture and the application of the method [] to randomize the low frequencies of the input. Statistics-based methods estimate a set of the sample's statistics, which are then enforced in the synthesized texture. The statistical characterization of texture images was initiated by Julesz. Julesz was the first to point out that texture images could be reliably organized, according to their N-th order statistics, into groups of textures that are preattentively indistinguishable by humans []. In [], Julesz demonstrated that many texture pairs sharing the same second-order statistics would not be discerned by human preattentive vision. This hypothesis constitutes the first Julesz axiom for texture perception. One consequence of this axiom is that two textures sharing the same Fourier modulus but with different phase should be perceptually equivalent. This motivates a class of algorithms (random phase methods) aiming at creating textures with a given secondorder statistic. An example of such algorithms is [] and its extension [] in which a texture is generated by randomizing the Fourier phase while maintaining the Fourier modulus. The random phase method in [] correctly synthesizes micro-texture images which adapt well to a Gaussian distribution, but it fails for more structured ones, as can be experimented in the executable paper []. It has then been showed that textures with the same second and even third order statistics were visually different [, ]. In [, ] Julesz proposed a second theory to explain texture preattentive discrimination theory by introducing the notion of textons (local conspicuous features). This new theory states that only the first order statistics of these textons are relevant for texture perception: images having the same texton densities (first ordre statistic) should not be discriminated. Relative to the texton theory, one axiom is that texture perception is invariant to random shifts of the textons [] leading to the stochastic dead leaves models [, , ].





List of Publications



Heeger and Bergen [] extended the Julesz approach to multiscale statistics. They characterized a texture sample by the histograms of its wavelet coefficients. By enforcing the same histograms on a white noise image, they obtained an exemplar based synthesis method. Yet this method only measures marginal statistics. It misses important correlations between pixels across scales and orientations. See the online execution of this method [] where some success but many failures are evident, like for []. Within a similar range of results De Bonet [] randomizes the initial texture image and preserves only a few statistics: the dependencies across scales of a multiresolution filter banch response. In [] the authors propose to synthesize a texture by taking randomly patches from the sample texture and placing them randomly in the output texture image. A blending step is applied across the overlapping blocks to avoid edge artifacts. The results achieved are similar to [, ]. Other methods are also based on statistics of wavelet coefficients or more involved multiscale image representations [, , ]. The Heeger-Bergen method was extended by Portilla and Simoncelli [] who proposed to evaluate on the sample some  crosscorrelations, autocorrelations and statistical moments of the wavelet coefficients. Enforcing the same statistics on synthetic images achieves striking results for a very wide range of texture examples. This method, which represents the state of the art for psychophysically and statistically founded algorithms is nevertheless computationally heavy, and its convergence is not guaranteed. Its results, though generally indiscernible from the original samples in a pre-attentive examination, often present blur and phantoms.

Patch-based methods constitute a totally different category of texture synthesis algorithms. The initial Efros and Leung [] method is based on Shannon's Markov random field model for the English language. In analogy with Shannon's algorithm for synthesizing sentences, the texture is constructed pixelwise. For each new pixel in the reconstructed image, a patch centered in the pixel is compared to all the patches of the input sample. The patches in the sample that are closer help predict the pixel value in the synthetic image. Several optimizations have been proposed to accelerate this algorithm. Among them Wei and Levoy [], who managed to fix the shape and size of the learning patch and Ashikhmin [] who proposed to extend existing patches whenever possible instead of searching in the entire sample texture. Yet, as already pointed out in the original paper [], an iterative procedure may fail by producing "garbage" when the neigborhood's size is too small. On the other hand it can lead to a trivial verbatim reproduction of big pieces of the sample when the neighborhood is too large. This can be experimented in []. Several extensions of [] have been proposed that manage to accelerate the procedure and resolve the "garbage" problem by stitching entire patches instead of pixels. In [] the authors propose to synthesize the image by quilting together patches that were taken from the input image among those who best match the patch under construction. A blending step was also added to overcome some edge artifacts. Efros and Freeman [] proposed an extension of the latter introducing for the quilting method a blending algorithm that computes a path with minimal contrast across overlapping patches, thus mitigating the transition effect from patch to patch. Kwatra et al. in [] extend [] by using a graphcut algorithm to define the edges of the patch to quilt in the synthesis image. Another extension of [] is proposed by Kwatra et al. in [] where to synthesize a texture image they improve the quality of the synthesis image sequentially by minimizing a patch-based energy function. For these methods the visual results are strikingly good, but they still retain the risk of verbatim copying large parts of the input sample. For practical applications this may result in the appearance of repeated patterns in the synthesized image.

An ideal exemplar-based texture synthesis algorithm should create a new texture image that is perceptually equivalent to the input sample, but this image should be as random as possible while maintaining the same statistics as the input. To the best of our knowledge, none of the methods



proposed in the literature are able to combine these advantages. On the one hand, statistics-based methods fail, in general, to be perceptually equivalent to the input, either because the statistics are not well enforced or because these statistics are not complete enough. On the other hand, patchbased methods are too close to a mere copy-paste. This effect can be observed in the synthesis maps of Figure .. In this chapter, we propose an algorithm that lies between statistics-based and patch-based algorithms. We gave up choosing a set of statistics from the input sample to enforce in the synthesized image. Instead, we chose to model each texture patch by a Gaussian distribution learned from its similar patches. Inspired by [], we maintain the idea of searching for patches to stitch together in the original sample. However, instead of using the exact patch taken in the input texture, we sample the stitched patch from its Gaussian model. Locally Gaussian patch model have been proved useful in image denoising []. Our approach permits to maintain the coherence between patches with respect to the input sample, while creating new patches that do not exist in the sample texture but are still perceptually equivalent to it. In this way, we manage to combine the positive aspects of statistics-based and patch-based methods, while overcoming some of their drawbacks.

This chapter is structured as follows. In Section ., we give an estimation method of the Gaussian model for each patch. An analysis of the model's variance is provided. In Section ., the locally Gaussian synthesis algorithm is described. Section . compares the results with state of the art texture-synthesis methods. The concluding section . also states the limitations of the method. = {u ((x, y) + (i, j)) , (i, j) ∈ [0, . . . , n -1] 2 }, (x, y) ∈ V u where V u = I M -n+1 ×I N -n+1 denotes the discrete domain of the valid patches in u. The patch p (x,y) u will be considered as a column vector of size dn 2 × 1 I m denotes the m × m identity matrix

. The patch Gaussian model

In this section we consider a first and simple model of the texture sample. We denote u : Ω → R the input texture image defined on the discrete domain Ω = I M × I N of size M × N and p (x,y) u a patch from u defined as a square block of size n × n among all possible overlapping patches in u. Each patch is defined by the position of its left top corner (x, y) ∈ Ω. We assume that the



underlying distribution of a given patch p (x,y) u in u follows a multivariate Gaussian distribution of mean µ (x,y) and covariance matrix Σ (x,y) . Hence to generate a new texture image w : Ω → R according to this model all the patches of u are replaced by samples of these multivariate Gaussian distributions.

Given the image u and fixing an overlap size o ∈ [0, 1] we will consider only the patches from u that are taken in a raster-scan order (left to right and top to bottom) and whose centers are separated by (1-o)n pixels. We call them the patches under construction. They can overlap in three possible way: horizontally, vertically and forming an L-shape overlap (see Figure .). Each one of these patches will be replaced by a new patch p(x,y) u sampled from its corresponding multivariate Gaussian ditribution of parameters µ (x,y) and Σ (x,y) .

The parameters µ (x,y) and Σ (x,y) of the Gaussian distribution of the patch p (x,y) u are estimated from the set of nearest patches U = {p

(x i ,y i ) u , i = 1, . . . , m}. Here p (x i ,y i ) u
are the m nearest patches to p (x,y) u in u according to the L 2 distance. The empirical statistics µ (x,y) and Σ (x,y) are then defined as

µ (x,y) = 1 m m i=1 p (x i ,y i ) u Σ (x,y) = m i=1 m j=1 p (x i ,y i ) u -µ (x,y) p (x j ,y j ) u -µ (x,y) t (m-1) . (.)
To simplify the notation, we denote by P the matrix whose columns are the normalized patches in vector form (p y) ), i = 1, . . . , m. Then we can reformulate the covariance matrix Σ (x,y) Σ (x,y) = 1 (m -1) P P t .

(x i ,y i ) u -µ (x,
Sampling step Sampling a vector p(x,y) u ∼ N (µ (x,y) , Σ (x,y) ) boils down to a linear combination of the m nearest patches p

(x i ,y i ) u such as p(x,y) u = 1 √ m -1 m i=1 a i (p (x i ,y i ) u -µ (x,y) ) + µ (x,y) = 1 √ m -1 P A + µ (x,y) (.)
where a i ∼ N (0, 1), i = 1, . . . , m are independent random variables and A is a column vector of size m × 1 whose elements are the standard normal independent random variables a i . We denote by N (µ, Σ) the multivariate normal distribution of mean µ and covariance matrix Σ. Indeed the sampled patch p(x,y) u follows the mutltivariate normal distribution N (µ (x,y) , Σ (x,y) ) where ,y) .

E(p (x,y) u ) = 1 √ m -1 m i=1 E(a i )(p (x i ,y i ) u -µ (x,y) ) + µ (x,y) = µ (x,y)  and E((p (x,y) u -µ (x,y) )(p (x,y) u -µ (x,y) ) t ) = 1 (m -1) m i=1 m j=1 E(a i a j )(p (x i ,y i ) u -µ (x,y) )(p (x j ,y j ) u -µ (x,y) ) t = 1 m -1 m i=1 (p (x i ,y i ) u -µ (x,y) )(p (x i ,y i ) u -µ (x,y) ) t = Σ (x
Blending step Once a patch p(x,y) u is sampled it is placed in image w at position (x, y) overlapping the part which has been synthesized previously. Across this overlap area the pixel values are blended using a convex combination. That is, at each position (i, j) in the overlap region two values are possible: the one from the previous synthesis g old and the one from the actual synthesis g new . We denote g the value assigned at position (i, j) defined as

g = (1 -γ(i, j))g old + γ(i, j)g new where γ = (1 -α)(1 -β) and α(x, y) = 1 -x-1 on-1 if 1 ≤ x ≤ on 0 if x > on and β = 1 -y-1 on-1 if 1 ≤ y ≤ on 0 if y > on .
The procedure is summarized in Algorithm . Compute U = {p y) , A ∼ N ( 0, I m ) {sampling step}

(x i ,y i ) u , i = 1, . . . , m} the set of m closest patches to p (x,y) u : Estimate µ (x,y) ← 1 m m i=1 p (x i ,y i ) u : p(x,y) u ← 1 √ m-1 (P -µ (x,y) )A + µ (x,



no longer represents the patch faithfully because the set U will contain outliers. Nevertheless, reasonable values for n and m ensure a correct reconstruction when replacing patches by others simulated from multivariate Gaussian distributions. We observed that when using 30×30 patches, it is reasonable to consider a neighborhood of 10 to 20 patches. Even though the generated texture is a random and new image perceptually equivalent to the input one, both are visually too similar. We would like to observe more randomness in the generated texture. A second inconvenient of this model is that the generated texture has necessarily the same size as the input image. For these reasons we propose a second model which does not fix the patches and has no constraint on the size of the generated texture. This is explained in the next section.

The model's variance

It is important to analyze how the variance of the Gaussian distribution varies with the patch size n and the number of nearest neighbours m used to estimate these distributions. We expect to have a variance that is not equal to zero since this would correspond to taking the best match in u and thus having no innovation at all. On the other hand a high variance would introduce too much variability and therefore end up with a blurry reconstruction of the texture.

The ranges for m and n yielding a correct estimate for the variance cannot be determined a priori



and strongly depend on the texture sample. To measure the variability of the Gaussian model we observe the mean standard deviation per pixel of a given patch p (x,y) u in the texture image which can be estimated by

σ = 1 dn 2 dn 2 i=1 Σ (x,y) (i, i) = 1 dn 2 dn 2 i=1 σ i
where d is the number of channels (d = 1 for grayscale images and d = 3 for color images).

In As expected for a fixed patch size n the mean standard deviation per pixel σ increases when using more patches to estimate the Gaussian distribution. For a fixed value of m the same behaviour is observed when increasing n. The first example, a micro texture, shows that the variance of the patch is not negligible even for m = 5. Gaussian models are well suited for this kind of texture. The second example, a periodic texture, shows that the patch mean standard deviation varies between 8 and 18 for values of m ∈ {5, 10, 30, 50}. These are reasonable values confirming that effectively the patches simulated have an acceptable degree of innovation.

Input 500 × 500



. The locally Gaussian texture synthesis algorithm

As mentioned in the previous section simply replacing the texture patches by Gaussian samples is not enough for synthesizing a texture. This model yields texture images that are visually too similar to the input sample and of fixed size. We aim at generating images that are perceptually equivalent to the input yet not visually identical. Inspired by the work of Efros and Freeman [] we propose an iterative algorithm that models the self-similarities of the input sample as multivariate Gaussian distributions.

The proposed algorithm consists in generating a new texture image w patch by patch in a raster-scan order given the input sample u. At each iteration the goal is to generate a new patch p (x,y) w that is partially defined on a region called the overlap region (see Figure .). The known part of the patch will define the set of patches from u that will infer its Gaussian model. The generated patch p We initialize the new texture image w placing a seed patch in its top left corner. For this we pick a random patch from the texture sample u, estimate its underlying Gaussian distribution and sample a new patch from it (the seed patch). In continuation, each synthesized patch is sampled from a Gaussian model learnt from the m patches with their left half most similar to the right half of the last synthesized texture patch. This new patch is then quilted on the current texture.

Patch stitching is an important step when adding the new patches at each iteration. In this chapter we propose to follow the procedure presented in []. In Chapters  and  we suggest other stitching alternatives. The method proposed in [] assigns new edges to the patches on the overlap area in order to quilt them along the path on which they best match. These edges are selected as curves minimizing the error between the new patch and part of the synthesized image across the overlap region. For a detailed description of this step please refer to Annex A... The synthesis procedure is summarized in Algorithm .

The proposed method yields a new texture image perceptually similar to the original sample and whose patches have been randomized. The result obtained is visually similar to the input sample, i.e. we could consider that the input texture and the new texture are two different sub images of one same bigger texture. Initializing the input with the seed patch and synthesizing the  Gaussian patches in a raster-scan order using only their overlap to find the nearest neighbours from the input yields this variability regarding the first model we proposed. A discussion on the underlying mathematical model The underlying random field of a given input texture is assumed to be a Gaussian random field. Indeed, as illustrated with the experiment in Figure ., a patch in the input is assumed to be a multivariate Gaussian vector. Furthermore, we assume the existence of a set of other patches following the same probability distribution, from which the parameters (mean and covariance matrix) are estimated. Thus the resulting reconstructed image is a sample of a Gaussian random field. Nevertheless, the same cannot be stated for the underlying random field of the synthesized textures when using Algorithm . These generated textures are not samples of a Gaussian random field. What can be affirmed is that the conditional distribution of the patches are Gaussian, meaning that if p 0 w , p 1 w , p 2 w . . . are the overlapping patches of the image taken in a raster scan order, then the conditional probability P(p n w |p n-1 w ) is a Gaussian multivariate distribution. However, the mean and the covariance matrix Initialize w placing a seed patch in its top-left corner (x = 1, y = 1). The image w is of size rM × rN , where M × N is the size of u.

: for x = 1 : n(1 -o) : rM do : for y = 1 : n(1 -o) : rN do : if (x > 1 or y > 1) then :
Find the m nearest patches p

(x i ,y i ) u
, i = 1, . . . , m in the input sample u that best agree with the patch under construction p (x,y) w along the overlap area

: Estimate µ (x,y) ← 1 m m i=1 p (x i ,y i ) u : Sample p(x,y) w ← 1 √ m-1 P A + µ ∼ N µ (x,y) , Σ (x,y) where A ∼ N ( 0, I m ) : Quilt p(x,y) w in w at position (x, y) {see Annex A..} : end if :
end for : end for : return w of this Gaussian distribution are computed from patches similar to p n-1 w in the input sample. Therefore the joint probability distribution of the patches is not necessarily Gaussian. Thus, while the algorithms are effective, they do not imply so far the existence of a complete texture model. Indeed, the input sample is characterized by a Gaussian distribution for each patch but also by spreads, defined as the spatial distribution of the patches belonging to the same Gaussian model. For some periodic textures the spread is deterministic (like in a chessboard for instance) but for general textures the spread itself is random. To model the whole texture as a random field a stochastic model for the spreads would therefore be necessary. For example, if the spread of the input sample is given and the patches of the generated texture are samples of conditional Gaussian distributions respecting exactly the same spread as the input's, then the underlying random field of the whole synthesized texture is Gaussian (as in the example of Figure .). Nevertheless, to synthesize textures this is not very interesting in terms of variability among the different sampled textures. Another meaningful example where the model is complete is the case of periodic textures where the spread is deterministic. For a fixed seed patch (the one that initializes the generated texture) the underlying random field of the generated textures is a periodic Gaussian random field. But this is no longer the case when the seed patch is random.

To summarize, we cannot claim that the proposed synthesis algorithm corresponds to a complete texture model, because it would require a stochastic model for the spread, namely the spatial interaction between patches obeying the same Gaussian model. This is left for future investigation, perhaps in the spirit of [].



. Experiments

A first set of synthesis results can be seen in Figure . to illustrate the influence of the parameters n and m for a given texture example. We vary the values of the number of nearest neighbours m used to estimate the parameters of the Gaussian distributions and of the side patch size n. The first observation is that the preferable patch size has a sample-dependent lower bound. If the patch size is too small the algorithm does not capture the variability of the sample. This is the case for all non parametric patch-based methods. Yet, the larger the patch size, the fewer the neighbours that we can use to build a faithful model for the patches. The second observation relates to the size of the overlap between patches. If the proportion of the overlap area with respect to the patch size is too small, the obtained model may no longer be a correct representation of the patch to simulate. Indeed, the patches are compared on the overlap area but the model is learnt on the whole patch. We observed that for an overlap of half the patch size the model is faithful enough to achieve good results. The number of nearest neighbours used to infer the Gaussian models is strongly dependent of the input texture. For periodic or pseudo-periodic textures this is not an issue, as long as the patch size is not too large, since there are many similar patches to infer a good patch model. Whenever strong and isolated configurations are visible then this becomes tricky and the value of m is limited. Moreover a unique value of m may not be adapted in particular if the texture is not stationary. Another limitation on the value of m is indeed the size of the input image.

Figure . compares several state of the art texture synthesis methods. In general we can observe that for statistics-based methods (three first columns) the quality of our algorithm's visual results is considerably superior. However, when the size of the patch n is too large compared to the texture's structure (last row in Figure .) blurring artifacts can appear in our results. On the other hand, like for the Efros-Freeman algorithm [] our method is able to preserve the strong structures available in the input sample. We obtain similar visual results while creating an image whose patches are new configurations. Of course one can notice that for isolated configurations even if the generated patches are pixelwise different to the input ones, visually they are very similar and verbatim copy effects are still visible. We believe that being able to work with smaller patches reduces significantly this drawback.

In Figure . we compare the capacity of innovation of our method to the work in []. To illustrate this we represent the synthesis results using position and synthesis maps as can be seen in Figure . which shows the input sample u, the associated position map (Pmap), the synthesized texture w and the associated synthesis map (Smap). Each pixel position x in the sample texture u is associated to a different color from a continuous colormap. The resulting image is the position map. The synthesis map associated to the synthesized texture w is obtained by mapping each of its positions x ′ to the color value of the corresponding position x of u. This position x corresponds to the central pixel of the nearest patch in u to the patch centered in x ′ in w. The synthesis map permits to identify the tendency of an algorithm to generate verbatim copies and to visualize from which regions of the input texture are sampled the patches. To compute these synthesis maps we used the PatchMatch algorithm [] an efficient method to approximate optimal correspondences. Patch match is an algorithm that finds approximate nearest neighbours matches between image patches. That is, given two input images A and B the patch match algorithm returns an offset map. This offset map indicates where, for each pixel x ′ i in B considering the patch centered in x ′ i , is approximatively the nearest patch in A centered in x j . We used the implementation of the patch match algorithm provided in [].



. Conclusion

In this chapter we presented a synthesis algorithm that generates a new texture image by stitching together patches sampled from multivariate Gaussian distributions in the input sample patch space. For each patch under construction in the output image the underlying Gaussian model is inferred from the set of similar patches (to the patch under construction) from the input sample. Hence the algorithm synthesizes a texture that is perceptually equivalent to the sample texture yet not composed of patches existing in the input texture. This method reduces some of the drawbacks of the statistics-based and the patch-based methods. The stitching procedure is a bit complex and will be replaced in Chapter . Like the Efros-Leung or the Efros-Freeman methods, the algorithm remains dependent on the choice of n and m, that may have to be adjusted for each texture sample. In our opinion the texture samples used in the literature are too small, particularly for macro-textures like the ones presented here and in most papers. Thus, our local Gaussian model is forced to use only  to  degrees of freedom because only some  to  patches are similar enough. This estimate should improve with larger texture samples. Our algorithm has low complexity, compared for instance with classic patch-based denoising algorithms [], []. Indeed, the patch covariance matrix needs not to be computed or inverted. The gaussian samples are simply obtained by linearly combining the nearest neighbours patches. Another alternative to reduce the dependency of the method to the patch size is to work in a multiscale approach. This extension is presented in the next chapter. Unlike the statistics-based algorithms, but like the other patchbased methods, our algorithm is not forced to respect the global statistics of the texture sample. This can be observed in the second row of Figure . where our result correctly reproduces an image of petals but fails to insert a correct proportion of the leaves in the synthetic image. In the next chapter we also deal with this issue.



Input

RPN

HB PS EF LG

 3 Conditional Gaussian model: a multiscale algorithm

Exemplar based texture synthesis is defined as the process of generating, from an input texture sample, new texture images that are perceptually equivalent to the input.

In the present chapter, we model texture self-similarity with conditional Gaussian distributions in the patch space in order to extend the use of stitching techniques. Then, a multiscale texture synthesis algorithm is introduced, where texture patches are modeled at each scale as spatially variable Gaussian vectors in the patch space.

The Gaussian distribution for each patch is inferred from the set of its nearest neighbours in the patch space obtained from the input sample. This approach is tested over several real and synthetic texture images and its results show the effectiveness of the proposed technique for a wide range of textures.

. Introduction

In this chapter, an extension and detailed description of a texture synthesis framework is proposed which combines a multiscale approach to the locally Gaussian texture model in the patch space. Each texture patch of the synthesized image is sampled from its Gaussian distribution estimated on a set of similar patches taken in the input sample as defined in Chapter . Modeling patches with conditional Gaussian distributions is an approach that is also used in image processing as for instance in the denoising algorithm of Buades et al. []. The first question that was brought up was how to stitch together the patches. In Chapter , the Efros and Freeman's [] was adopted, where every new patch overlaps the previously synthesized one, and the overlapped parts are blended.

In this chapter we consider conditional Gaussian models constrained to respect the values of its overlapping zone. This solution shows satisfying results for periodic or pseudo-periodic textures but fails for more complex ones. The second and central question brought up by non paramatric patch based methods is: how to handle the strong dependency of the method on the patch size in particular when dealing with macro-textures. Macro-textures show information at different scales that cannot be captured with a unique patch size. To this aim, a multiscale approach is considered. A full description and discussion is provided here.

In this chapter, we aim at extending the Gaussian model proposed in Chapter  to take into account the stitching step. We first recall briefly the principle of non-parametric patch-based methods and review the most recent contributions. As we have seen in Chapter , non-parametric patch-based methods were initialized by Efros and Leung [] who extended to images Shannon's Markov random field model initially devised to simulate text. In analogy with Shannon's algo- rithm for synthesizing sentences, the texture is constructed pixelwise. For each new pixel in the reconstructed image, a patch centered at the pixel is compared to all patches with the same size in the input sample. The nearest matches in the input help predict the pixel value in the reconstructed image. Several optimizations have been proposed to extend and accelerate this algorithm. Among them, Wei and Levoy [], managed to fix the shape and the size of the learning patch and Ashikmin [] proposed to extend existing patches whenever possible instead of searching in the entire sample texture. As we mentioned, these pixelwise algorithms are not always satisfactory. They are known to grow "garbage" when the compared patches are too small or when the input texture is not stationary, or may lead to verbatim copies of significant parts of the input sample for large patches as can be verified in []. This explains why more recent methods stitch together entire patches instead of performing a synthesis pixel by pixel. The question then is how to blend a new patch in the existing texture. In [] this is done by a smooth transition. Efros and Freeman [] refined this process by stitching each new patch along a stitching path with minimal contrast between the new patch and the texture in their overlapping zone. Kwatra et al. in [] extended the stitching procedure of [] by a graph cut approach redefining the edges of the patches. In [] the authors proposed to synthesize a texture image by sequentially improving the quality of a synthesis by minimizing a patch-based energy function. In the same spirit as [], where texture optimization is performed, the authors in [] proposed to synthesize textures in a multiscale framework using the coordinate maps of the sample texture at different scales. They introduced spatial randomness by applying a jitter function to the coordinates at each level combined to a correction step inspired by []. One of the key strengths of the method is that it is a parallel synthesis algorithm which makes it extremely fast. These non-parametric patch-based approaches often present satisfactory visual results. In particular they have the ability to reproduce highly structured textures. However, the risk remains of copying even several times verbatim large parts of the input sample. Furthermore, a fidelity to the global statistics of the initial sample is not guaranteed, in particular when the texture sample is not stationary. We refer to [] for an extensive overview of the different non-parametric patch-based methods.

Recently methods such as [, ] combine patch-based and statistics-based methods to overcome the drawbacks mentioned previously. In [], the author proposes to use a patch-based approach where all the patches of the synthesized image are created from a sparse dictionary learnt on the input sample. In [], Tartavel et al. extend the work of [] by minimizing an energy that involves a sparse dictionary of patches combined to constraints on the Fourier spectrum of the input sample. Still more recently, Gatys et al. [] have introduced the neural network methodology into the field. Extending the parametric approach, they characterize a texture by the cross-correlations of kernels learnt from a convolutional neural network dedicated to shape recognition. This algorithm emulates complex textures containing large objects for which the Portilla and Simoncelli algorithm is not satisfactory. Contrarily to the Portilla and Simoncelli approach, no Occam's razor was applied to reduce the number of texture parameters. The number of kernel correlations involved is perhaps exceedingly large: it appears that parts of the input are being recombined in the output. Nevertheless this learning approach is definitely promising.

The rest of this chapter is structured as follows. In Section . the local Gaussian (LG) model is described. An analysis of the model's variance is provided as well as the description of a first synthesis method. In Section . two new patch models are introduced: the conditional local Gaussian model (CLG) and the regularized conditional local Gaussian model (RCLG). Section . presents the multiscale approach (MSLG). Section . shows several experiments: a comparison of the three patch models proposed, a comparison to the results of state of the art texture synthesis methods and the influence of the parameters involved in the multiscale texture synthesis algorithm.

-1) u k Ω k → R: zoom out of u of a factor 2 k , defined on the discrete domain Ω k = I 2 -k M × I 2 -k N of size 2 -k M × 2 -k N for k = 1 . . . K -1 w k Ω k r → R: synthesized texture at scale k, defined on the discrete domain Ω k r = I r2 -k M × I r2 -k N of size r2 -k M × r2 -k N for k = 0 . . . K -1 v k Ω k r → R: zoom in of w k+1 of factor 2. It is the initialization of the low resolution of the synthesized image w k for k = 0 . . . K -2 p (x,y) u square patch of size n × n from an image u of size M × N at position (x, y), p (x,y) u = {u ((x, y) + (i, j)) , (i, j) ∈ [0, . . . , n -1] 2 }, (x, y) ∈ V u where V u = I M -n+1 ×I N -n+1
denotes the discrete domain of the valid patches in u. The patch p (x,y) u will be considered as a column vector of size dn 2 × 1

G σ

Gaussian kernel of mean zero and standard deviation σ, G σ (x, y) = 1 2πσ 2 e -x 2 +y 2 2σ 2

L u k Ω k → R: low resolution of image u k , L u k = u k * G σ , k = 0 . . . K -2 H u k Ω k → R: high resolution of image u k , H u k = u k -u k * G σ , k = 0 . . . K -2 L w k Ω k r → R: low resolution of image w k , L w k = w k * G σ , k = 0 . . . K -2 H w k Ω k r → R: high resolution of image w k , H w k = w k -w k * G σ , k = 0 . . . K -2 I m
denotes the m × m identity matrix

. Gaussian patches

For all notations used in this section we refer to Table . for a detailed definition.

.. The assumption of a Gaussian patch model

We recall in this section the locally Gaussian model suggested in Chapter . For a given texture image u, the underlying distribution of every patch p (x,y) u is modeled as a multivariate Gaussian distribution of mean µ (x,y) and covariance matrix Σ (x,y) .

These parameters are estimated from the set of nearest patches U = {p

(x i ,y i ) u , i = 1, . . . , m} where p (x i ,y i ) u
are the m nearest patches to p (x,y) u in u according to the L 2 distance. The empirical  statistics µ (x,y) and Σ (x,y) are then defined as

µ (x,y) = 1 m m i=1 p (x i ,y i ) u Σ (x,y) = 1 (m-1) P P t . (.)
We denote by P the matrix whose columns are the normalized patches in vector form (p

(x i ,y i ) u - µ (x,y) ), i = 1, . . . , m.
The sampled vector p(x,y) u is thus defined as a linear combination of the m nearest patches p

(x i ,y i ) u p(x,y) u = 1 √ m -1 P A + µ (x,y) . (.)
Here A is a vector whose elements are standar normal independent random variables A ∼ N ( 0, I m ). Thus, the sampled vector has the correct Gaussian distribution N (µ (x,y) , Σ (x,y) ).

In Chapter  we measured the variability of the Gaussian models with respect to the patch size n and the number of paches m. For this the mean standard deviation per pixel of a given patch in the texture image was estimated as

σ = 1 dn 2 dn 2 i=1 Σ (x,y) (i, i) = 1 dn 2 dn 2 i=1 σ i
where d is the number of channels(d= for grayscale images and d = 3 fop color images).We observed that indeed the models considered have reasonable variances confirming that effectively the patches simulated have an acceptable degree of innovation. From now on to simplify the notations when referring to µ (x,y) and Σ (x,y) they will be denoted by µ and Σ respectively.

.. Texture synthesis algorithm

We remind the locally Gaussian texture synthesis algorithm (LG) presented in Chapter  with a different notation that will be useful for the following section. Given the input texture u, an output image w is synthesized sequentially patch by patch in a raster-scan order (left to right, top to bottom). Each new patch added to w overlaps part of the previously synthesized patch as can be seen in Figure .. Depending on the stage of the synthesis three different cases of overlap can be observed: vertical (first row of raster-scan), horizontal (first column of raster-scan) and L-shape (everywhere else). This is also shown in Figure .. Each patch is simulated following a multivariate Gaussian distribution of mean µ and covariance matrix Σ as defined in (.). When quilting the patch in w the same stitching step is applied as in Efros and Freeman's work [] where the authors propose to compute an optimal boundary between the new patch and the previously synthesized one along their overlap region. This is done thanks to a linear programming optimization.

To define the set U = {p

(x i ,y i ) u , i = 1, . . . , m} of nearest patches to p (x,y) w , let us consider p (x,y) w
as the patch being currently synthesized and taken as a column vector of size dn 2 × 1. The patch p (x,y) w will overlap part of the previous synthesis. To synthesize p (x,y) w

, we decompose it as

p (x,y) w = S t Sp (x,y) w + R t Rp (x,y) w . (.)
Here, S : R dn 2 → R dn 2 -k and R : R dn 2 → R k are projection operators such that Rp . This decomposition will be useful in the following section.

The patches p

(x i ,y i ) u
used to learn the parameters of the multivariate Gaussian distribution (.) are the m nearest neighbours in u to the current patch p (x,y) w , for the L 2 distance restricted to the overlap area, given by Rp In Figure . the synthesis results using the Gaussian sampling are compared to those of the quilting method [] to illustrate that the Gaussian sampling achieves visual results that are comparable to those in [] while providing a local parametric model, as shown in Chapter . With the Gaussian sampling some blur is introduced in the synthesis results but the effects of verbatim copies and garbage growing are reduced as can be observed in the second example in Figure ..

Regarding the uderlying mathematical model of the synthesized image, as discussed in Chapter , we cannot claim that the proposed synthesis algorithm corresponds to a complete texture model, because it would require a stochastic model for the spread, namely the spatial interaction between patches obeying the same Gaussian model.

. Conditional Gaussian models

The purpose of this section is to analyze how to replace the stitching step by a conditional patch model constrained to the overlap zone. The natural idea is to condition the new Gaussian samples to the values observed in the overlap zone, thus allowing to maintain the same or slightly modifying the overlap pixels. We observed that conditional Gaussian models used in [] as an inpainting method for microtextures achieved good results.

This requires modeling the transition effect between patches. Each new patch will be estimated as a Gaussian vector conditioned to the pixel values of its corresponding overlap region. In this way the simulated patch would naturally "agree" with w in the overlap area, thus avoiding a stitching procedure.

Such models will be considered in the following sections.



Input

LG method []

.. A first conditional Gaussian model

In this first patch model the idea is to model sample patches obeying a Gaussian distribution that exactly match with their overlap area pixels, avoiding the unwanted discontinuities. We shall call this model Conditional Locally Gaussian (CLG). Each patch p (x,y) w is taken as a column vector of size dn 2 × 1 and can be partitioned as

p (x,y) w = S t Sp (x,y) w + R t Rp (x,y) w .
We assume throughout that the vector p (x,y) w follows a Gaussian distribution of mean µ and covariance matrix Σ. These parameters can be expressed as follows :

µ =S t Sµ + R t Rµ = S t µ 1 + R t µ 2 and Σ =S t SΣS t S + S t SΣR t R + R t RΣS t S + R t RΣR t R =S t Σ 1,1 S + S t Σ 1,2 R + R t Σ 2,1 R + R t Σ 2,2 R.
The problem can be formulated as finding the "best sample" p = (p 1 , p2 ) conditioned to the overlap values q 0 that are known, p1 = arg max

p 1 P µ,Σ (P 1 = p 1 | P 2 = q 0 ). (.)



Here P 1 corresponds to the unknown values of p (x,y) w and P 2 to the values of p (x,y) w on the overlap area, i.e. P 1 = Sp (x,y) w and P 2 = Rp (x,y) w . Using classic results on conditional multivariate Gaussian distributions [], the distribution of P 1 conditioned to P 2 = q 0 is a multivariate Gaussian distribution of parameters μ and Σ where

μ = µ 1 + (SΣR t )(RΣR t ) -1 (q 0 -µ 2 ) and Σ = (SΣS t ) -SΣR t (RΣR t ) -1 RΣS t .
Since the most probable sample of a multivariate Gaussian distribution is its mean, the solution to (.) is

pw = R t q 0 + S t (µ 1 + (SΣR t )(RΣR t ) -1 (q 0 -µ 2 )). (.)
This first result has the advantage that since the mean of the conditional distribution is taken, the risk of unwanted sharp transitions between patches for some types of textures will be attenuated. Another observation is that since this new model might not have very probable samples, working with the most probable one improves the visual aspect of the sampled patch. The negative aspects of taking the most probable sample is that on the one hand the result, for a given initialization, is deterministic. On the other hand the synthesis result loses details by displaying the mean of the conditional distribution.

To overcome both mentioned drawbacks we sample a patch from the underlying Gaussian distribution, that is,

p = R t q 0 + S t p1 , p1 ∼ N (μ, Σ). (.)
Yet equations (.) and (.) show that these solutions do not make sense if (RΣR t ) is not invertible. This is unfortunately frequent, as the number of neighbours m used to build the Gaussian distribution is often very small compared to the dimension of the vectors we aim to model. Therefore the learnt Gaussian models are strongly degenerated. The fact that Σ is not invertible implies that the Gaussian vectors p (x,y) w ∼ N (µ, Σ) live in a subspace of R dn 2 . However it does not say that there is no solution to (.). Indeed there are two possible cases:

. The subspace of Gaussian vectors (N (µ, Σ)) intersects the subspace of vectors R t q 0 + S t p 1 .

Thus the Gaussian distribution N (μ, Σ) exists and there is a solution to our problem.

. The subspace of Gaussian vectors (N (µ, Σ)) does not intersect the subspace of Gaussian vectors R t q 0 + S t p 1 . Therefore there is no solution.

To avoid the second case (having no solution) a small perturbation is added to the Gaussian distribution learnt for p (x,y) w as follows

p (x,y) w ∼ N (µ, Σ + σ 2 I dn 2 ),
where σ 2 is a real positive number and I dn 2 is the identity matrix of size dn 2 × dn 2 ensuring that the problem is well conditioned. In that way the Gaussian vectors p ∼ N (µ, Σ) intersects the set of Gaussian vectors R t q 0 + S t p 1 , using the new distribution N (µ, Γ) will slightly modify the solutions in (.) and (.) when using small values of σ 2 . Thus it is enough to take a low value for σ 2 and the solutions obtained in both cases (N (µ, Σ) and N (µ, Γ)) will be very close to each other.

Adding this perturbation to the multivariate Gaussian model has the drawback of increasing the computational cost of sampling a Gaussian vector. For the model in Section .. a Gaussian vector was sampled by a simple linear combination of m patches. There was no need in estimating the covariance matrix of the Gaussian distribution. The model presented in this section requires the computation of the Gaussian parameters, in particular the covariance matrix Γ = (SΓS t ) -SΓR t (RΓR t ) -1 RΓS t . The covariance matrix Γ is symmetric and definite positive and then admits a Cholesky factorization, i.e. there is a unique lower triangular matrix C such as Γ = CC t . Thus the vector s = Ct + μ yields a sample of the desired Gaussian distribution N (μ, Γ) with t ∼ N ( 0, I dn 2 ). The complexity of computing the Gaussian samples with the cholesky decomposition is O(n 3 ) whereas the complexity of the first model is O(nm).

This synthesis method is summarized in Algorithm  and Algorithm .

Algorithm : CLG sampling Input: set of nearest pathces U , overlap values q 0 , projection operators R and S Output:

Gaussian sample F : µ ← 1 m m i=1 p (x i ,y i ) u
, mean of the Gaussian model : Σ ← 1 (m-1) (P P t ), covariance matrix of the Gaussian model

: Γ ← Σ + θ 2 I n 2 , covariance matrix of full rank : μ ← Sµ + (SΓR t )(RΓR t ) -1 (q 0 -Rµ), conditional Gaussian mean : Γ ← (SΓS t ) -SΓR t (RΓR t ) -1 RΓS t , conditional covariance matrix : C ← cholesky decomposition of matrix Γ : t ← Gaussian sample N ( 0, I n 2 ) : s ← Ct + μ Gaussian vector following N (μ, Γ) : return F

.. The regularized conditional Gaussian model

In the previous section the patches' statistical models were conditioned to the exact values of the synthesized pixels across their overlap area. This is too restrictive for some types of textures and is at risk of creating unlikely samples. Instead of forcing each patch p (x,y) w to take the exact same values on the previously synthesized part, it is therefore natural to allow the patch to vary slightly on the overlap area. This variation is rendered necessary by the scarcity of patch samples in a small texture sample. We shall call this model Regularized Conditional Locally Gaussian (RCLG). Consider the same patch model N (µ, Σ), but let us now allow the overlap components P 2 = Rp (x,y) w to take values P 2 = q 0 + b, where b ∼ N (0, θ 2 I k ) and where θ is the variation allowed for the overlap values. The joint distribution of (P 1 , P 2 ) is defined by:

for y = 1 : n -o : rN do : if (x > 1 or y > 1) then : q 0 ← known overlap values : U ← set of m nearest neighbours to p (x,y) w : p(x,y) w ← CLG sampling(U , q 0 ) : Quilt p(x,y)
w in w at position (x, y)

: end if :
end for : end for : return w

P(p 1 , p 2 ) =P µ,Σ (p 1 |p 2 )P q 0 ,θ 2 I k (p 2 ) = 1 Z e -1 2 ((p 1 ,p 2 )-(μ,q 0 )) t ∆ -1 ((p 1 ,p 2 )-(μ,q 0 )) , (.)
where Z = 2π det(∆) and ∆ = Σ 0 0 θ 2 I k .

Proceeding as in the previous model (CLG) the most probable sample p = arg max

(p 1 ,p 2 ) P µ,Σ (p 1 |p 2 )P 0,θ 2 I k (p 2 -q 0 )
is exactly the same as in (.) and has the same drawbacks. On the other hand sampling from (.) allows to relax the overlap constraint.

Once again, to guarantee the existence of the solution, the Gaussian distribution of p (x,y) w is slightly modified as done in the CLG model. A small perturbation is added to the covariance matrix Σ and the problem is then well conditioned. Thus, as in the CLG model, the computational cost of sampling from this new Gaussian distribution is higher compared to the LG sampling as mentioned for the CLG model since once again the parameters of the RCLG are estimated. The sampling process is similar to the one described in Algorithm . The difference is that the known values of the overlap are slightly modified to relax the constraint. Apart from that the procedure is the same where in one case the overlap values are exactly those from the previous synthesis (CLG) and in the second case they are slightly modified as q 0 + b where b ∼ N ( 0, θ 2 I k ).

. A multiscale generalization

For all notations used in this section we refer to Table . for a detailed definition.



Macro textures present details at different scales: a coarse one that contains the global structure of the texture and finer ones containing the details. Small patch sizes may capture the finer details of the input but the resulting texture will lack global coherence. On the other hand using large patches will maintain better the global structures on the risk of a "copy-paste" effect. Furthermore with large patches it becomes impossible to model the patch variability due to the curse of dimensionality, in other terms the lack of sufficient samples. This is for example apparent in the examples of Chapter  where modeling patches as multivariate Gaussian vectors leads to a slightly blurry texture.

Multiscale approaches permit to contemplate several patch sizes within one synthesis, i.e. to capture the different levels of details. If we fix the patch size to be n × n and use K scales this is similar to using K different patch sizes going from 2 K-1 n × 2 K-1 n to n × n for the coarsest to finer details within one synthesis.

In this section the potential of a multiscale approach is illustrated by improving the method described in Algorithm . Let us first introduce some notations. The input texture sample is denoted by u and u k , k = 1, . . . , K -1 are the zoomed out versions of u by a factor 2 k , k = 1, . . . , K -1. The synthesis result at each scale is denoted by w k , k = 1, . . . , K -1 and w is the synthesis result returned by the multiscale algorithm. An additional image needed at each scale is the low resolution of the result w k that is denoted by v k and its the result of zooming in w k+1 . The multiscale approach can be summarized in a few sentences. The method begins by a synthesis at the coarsest scale (k = K -1) using the local Gaussian method in Algorithm  where the quilting step is replaced by a simple average of the overlapping patches. For the remaining scales (k = K -2, . . . , 0) a synthesis is performed by using the result of the previous scale (k+1) and the input of corresponding resolution. At each scale the synthesis is done patch by patch in a raster-scan order. Each new patch, added to the synthesized image, overlaps part of the previously synthesized patch and is the combination of a low resolution patch and a high resolution one sampled from a multivariate Gaussian distribution. The Gaussian distribution of the high frequencies of a given patch is estimated from the high frequencies of its m nearest neighbours in the corresponding scale input image. The synthesis result of the finer scale is the desired output image.

We shall call this method Multiscale Locally Gaussian (MSLG). In the following the different parts of the method described in Algorithm  are detailed.

Zoom-out

The zoom-out operation is based on Gaussian filtering. It is performed as a smooth frequency cut-off followed by a sub-sampling of factor two. These operations are detailed below in (.), (.) and (.).

The smooth frequency cut-off is performed with the Gaussian kernel

G σ (x, y) = 1 σ 2 2π e -x 2 +y 2 2σ 2 , (.)
where we chose σ = 1.4. Indeed as proved in [] if we consider a Gaussian kernel of standard deviation σ = (0.8α) 2 -0.8 2 to smooth the image and then sub-sample it by a factor α we can assume that no aliasing is generated.

The sub-sampling operation by a factor two applied to an image u :

I M × I N → R is defined by S 2 : I M × I N → I M/2 × I N/2 , S 2 (u)(i, j) = u(2i, 2j), (i, j) ∈ I M/2 × I N/2 . (.)



The zoom-out of an image u of factor two is then

Z out 2 (u) = S 2 (u * G σ ). (.)
The images obtained after applying this zoom-out operation have no aliasing and they are blurry enough to avoid ringing artifacts when applying the zero padding zoom-in.

Zoom-in For an image u :

I M × I N → R, v = Z in 2 (u) : I 2M × I 2N → R
is the zoom-in by a factor two. This operation is performed by a zero padding of û : Î2M × Î2N → C where û denotes the discrete Fourier transform of u and Îc denotes the discrete interval [-c/2, . . . , c/2 -1] for c even and [-(c -1)/2, . . . , (c -1)/2] for c odd. This is done as

v(ξ 1 , ξ 2 ) = û(ξ 1 , ξ 2 ) if |ξ 1 | ≤ M -1 2 , |ξ 2 | ≤ N -1 2 0 else (.)
where v is the discrete Fourier transform of v and ⌊x⌋ denotes the integer part of x.

The relation between the zoom in and zoom out operators

Z in 2 (Z out 2 (u)) = u is valid when the image u verifies û(ξ 1 , ξ 2 ) = 0, ∀|ξ 1 | > M -1 2 , |ξ 2 | > N -1 2 .
We could have chosen other interpolation techniques as for example a spline interpolation. But a zero padding is well suited due to the nature of the zoomed out images which are blurry enough to avoid any ringing artifacts.

Distance between patches

To estimate the parameters of the Gaussian distribution of the patch being processed, denoted by p

(x ′ ,y ′ ) w k , the set U of m nearest patches in u k to p (x ′ ,y ′ ) w k
is considered. These patches are those minimizing the distance to p (x ′ ,y ′ ) w k defined in (.) for k = K -1 and in (.) for the remaining scales k = K -2, . . . , 0.

The size of patch overlap is fixed to half the patch size n/2. As mentioned in Section . there are three overlap cases: vertical (V.O.), horizontal (H.O.) and L-shape (L.O.). Here they are denoted as

Op (x,y) u = {u ((x, y) + (i, j)) , (i, j) ∈ O} where O =            [0, . . . , n -1] × 0, . . . , n 2 -1 if V.O. 0, . . . , n 2 -1 × [0, . . . , n -1] if H.O. 0, . . . , n 2 -1 × [0, . . . , n -1] ∪ n 2 , . . . , n -1 × 0, . . . , n 2 -1 if L.O. When k = K -1, the m nearest neighbours in u K-1 to the current patch p (x ′ ,y ′ )
w K-1 are those minimizing the L 2 distance restricted to the overlap area:

 d(Op (x,y) u K-1 , Op (x ′ ,y ′ ) w K-1 ) 2 = 1 |O| (i,j)∈O (u K-1 (x + i, y + j) -w K-1 (x ′ + i, y ′ + j)) 2 . (.)
When k = K-2, . . . , 0, the nearest neighbours in u k to the patch p

(x ′ ,y ′ ) w k
are those minimizing a distance (.) similar to (.) with an additional term taking into account the low resolution v k (the synthesis result of the previous scale k + 1). In (.) L u k denotes the low resolution of the image

u k , L u k = u k * G σ , k = 0, . . . , K -2. It is important to notice that when comparing Op (x,y) u k and Op (x ′ ,y ′ ) w k
the low and the high resolution must be considered jointly, they are not independent.

d(p (x,y) u k , p (x ′ ,y ′ ) w k ) 2 = 1 |O| (i,j)∈O (u k (x + i, y + j) -w k (x ′ + i, y ′ + j)) 2 + 1 n 2 n-1 i,j=0 (L u k (x + i, y + j) -v k (x ′ + i, y ′ + j)) 2 (.)

Blending process

The blending process consists in simply averaging the values across the overlap area as in (.). This step is applied only for the synthesis of scale k = K -1.

w k (x + i, y + j) = 1 2 (p (x,y) w K-1 (i, j) + p (x,y) w K-1 (i, j)) if (i, j) ∈ O p(x,y) w K-1 (i, j) if (i, j) ∈ I 2 n -O (.)
Synthesizing patches at scales k = K -2, . . . , 0 At each scale k a patch p (x,y)

w k is synthesized as the combination of a low resolution patch with a high resolution one. It can be decomposed as

p (x,y) w k = p (x,y) w k * Gσ + (p (x,y) w k -p (x,y) w k * Gσ ) = p (x,y) v k + (p (x,y) w k -p (x,y) v k ) = p (x,y) Lw k + p (x,y) Hw k .
Here L w k denotes the low resolution image of w k defined as

L w k = w k * G σ , k = 0 . . . K -2 and H w k denotes the high resolution image of w k defined as H w k = w k -w k * G σ , k = 0 . . . K -2.
The set U defines the Gaussian distribution of p

(x,y) Lw k ∼ N (µ L , Σ L ) and p (x,y) Hw k ∼ N (µ H , Σ H ) and
therefore the distribution of the patch p (x,y) v k . In this way the correlations between high and low resolution pixels are respected, using the low resolution synthesis v k as initialization.

w k ∼ N (µ, Σ) where µ = µ H + µ L and  Σ = Σ L + Σ H + cov p (x,y) Lw k , p (x,y) Hw k + cov p (x,y) Lw k , p (x,y) Hw k t . Instead of sampling p (x,y) Lw k from its Gaussian distribution, p (x,y) v k ∼ N (µ L , Σ L ) is kept
For all scales beside the coarsest one the synthesis is done by adding the high frequencies of the corresponding scale on the the low resolution basis image. It is the important to achieve a correct basis image in the coarsest scale on which the high frequencies will be added. The texture synthesis method is summarized in Algorithm .

Algorithm : Multiscale texture syntesis algorithm (MSLG)

Input: input texture sample u, side patch size n, number of nearest neighbours m, number of scales K, ratio (output size)/(input size) r Output: synthesized texture w

: Define u k ← Z out 2 (u k-1 ), k = 1 . . . K -1 : Define L u k ← u k * G σ , k = 0 . . . K -2 : Synthesize w K-1 ← LG(u K-1 , n, m, r) {see Algorithm } : for k = K -2 : 0 do : v k ← Z in 2 (w k+1 ) :
Initialize w k with zeros of same size as v k :

for x = 1 : n/2 : (r2 -k M -n + 1) do :
for y = 1 : n/2 : (r2

-k N -n + 1) do :
Compute d(p

(x ′ ,y ′ ) u k , p (x,y) w k ), for all (x ′ , y ′ ) in V u k , where V u k denotes the discrete domain of the valid patches in u k . : U ← {p (x i ,y i ) u k , i = 1, . . . , m} set of m nearest patches in u k that minimize d(p (x ′ ,y ′ ) u k , p (x,y) w k ) : H ← {p (x,y) u k -p (x,y) u k * G σ , ∀p (x,y) u k ∈ U}, high frequency of the corresponding patches in U :
Learn (µ, Σ) the parameters of the multivariate Gaussian distribution on the patches of H :

Sample p ∼ N (µ, Σ) : p(x,y) w k ← p (x,y) v k + p : w k ((x, y) + (i, j)) ← p(x,y) w k for (i, j) in I 2 n :
end for :

end for : end for : return w



. Experiments

In this section, texture synthesis results are shown using the algorithm described in Algorithm . In In general the results shown in Figure . are satisfying for a wide range of textures. The global structures are reproduced by the multiscale approach, while the local structures are maintained by the patch based approach. Using a Gaussian patch model allows to create new patches that do not exist in the input example while maintaining satisfying visual results. Based on the analysis of the patch's variance illustrated in Section . of Chapter  this guarantees that indeed the simulated patches are sufficiently different from the ones in the input sample. The examples of the last two rows illustrate some failure cases of the method. The main failure cause is the size of the input texture. It is obvious that the input must be large enough to provide us with a sufficient number of patch samples to estimate their distribution. If that is not the case, even though the patches are pixel-wise different of the input ones the visual aspect may remain too similar and cause a "copypaste" effect when m is small enough. Furthermore, as in other non parametric methods,"garbage", namely the excessive use of a subset of patches in the input image is not fully avoided. This effect is nevertheless mitigated by the multiscale approach.

.. Model comparison

The results in Figure . show the effects of avoiding the use of the blending step in Algorithm . The three models LG, CLG and RCLG were tested on several types of textures. For the three of them the quilting step was omitted. This was done to achieve a better comparison of the capacity of respecting the overlap by modeling or not the restriction of the patches to the overlap values. The results are compared for the three models applied in the sampling mode. The results of the same models in the best sample patch mode are less interesting and are not shown here. The general conclusion of the results in Figure . is that the conditional patch models achieve a better transition between patches on their overlap region, as expected. Nevertheless the results for CLG and RCLG both loose some resolution compared to the LG results. One can observe in Figure . that blur appears progressively in the raster-scan order of the synthesis algorithm. Indeed, when moving forward in the synthesis the m nearest patches start being too different from each other (due to the strong restriction of keeping the values of the overlap) and the result is a much blurrier patch. The first five texture examples in Figure . show better stitching results for CLG and RCLG than for to LG. Looking carefully at the synthesis results of the LG model the edges between patches are more noticeable than in the other two models where almost no transition effect can be seen. Nevertheless the loss of resolution caused by the Gaussian distributions is increased in both conditional models, mostly for CLG. In the last row example the results are less convincing. There are not enough reliable samples to estimate a correct conditional model. The visual results of CLG and RCLG get more and more degraded in the order of the raster scan. This is not surprising since the sample patches respecting the overlap values become increasingly unlikely.

From the experiments in



for complex textures and therefore the LG model is better to model locally the texture input. It generates less blur in the results and has a significant smaller computational cost. In the rest of the experiments only the LG model in considered. Nevertheless it could be interesting to test the use of conditional models in the multiscale version at the coarsest scale. In Figure . the LG model is compared to the MSLG model. One can observe the strength of using the multiscale approach in terms of reproducing global arrangements that are not kept in the LG model for a fixed patch size. This allows to obtain satisfying results for small values of n and therefore gives more freedom to choose the number of nearest patches m. One can also conclude from the experiments in Figure . that the multiscale approach generates blurrier textures than the LG model.

.. Comparison to other texture synthesis methods

In this section the results of Algorithm  are compared to other synthesis methods such as [, , ]. In general, the results obtained with the multiscale locally Gaussian method are visually comparable to the non-parametric patch based method of Efros and Freeman [], with the advantage that now, the patches are being sampled from their Gaussian model and therefore are different from their original patches. In Figure ., the first column shows the result of the proposed method. A noticeable drawback of the method is the loss of resolution caused by the use of Gaussian distributions. The proposed method was therefore combined with the Portilla and Simoncelli's algorithm [] as a first and simple approach to overcome this resolution loss. To combine both methods we first synthesize the texture example with MSLG and use this synthesis result as the initialization of the Portilla and Simoncelli algorithm instead of a "noise image". The result can be seen in the second column of Figure .. The result of combining both algorithms is very satisfying. On the one hand the local and global structures are kept due to use of the patch based and multiscale method. On the other hand [] allows to respect the global statistics of the input or at least be quite close to them. Comparing columns one and two to the third one shows how the combination of both methods improves the results of each method used separately. Of course this solution is limited in particular when the size of the synthesized image is too large. The fourth column shows the results of the Tartavel et al. method []. It is interesting to compare our results to this method since both approaches are multiscale, patch based, and create systematically new patches. The results for organized highly structured textures are comparable for both cases. Nevertheless for more complex textures like the flower example and the last two rows one notices a lack of sharpness when recreating the salient objects of the input with the method in []. Finally the last column shows results of []. One can observe that for MSLG the visual results are in general comparable to those of [] while providing a local parametric model. The results of [] are obviously excellent. But once again, in this kind of method, the algorithm ends up copying very large parts of the input. To illustrate this we represent our synthesis results and the ones obtained with [] using the position and synthesis maps in section . of Chapter . The synthesis map permits to identify the tendency of an algorithm to generate verbatim copies and to visualize from which regions of the input texture are sampled the patches. To compute these synthesis maps we used the PatchMatch algorithm []. In Figure . one can observe that in general the synthesis maps associated to the results of MSLG are more "noisy" than those associated to the results of []. Also in the synthesis maps associated to [] larger continuous zones are identified. This corresponds to the verbatim copies produced by the method. It is important to notice, that in both cases, for some texture examples in the synthesis maps there are some dominant colors represented in the synthesis maps. That reflects the discussion of Section . where we said that the spread of the input sample is not being



Input

LG model CLG model RCLG model



Input

LG model MSLG model



respected. For example, in the fifth row of Figure ., the input sample is not stationary (there is a change of luminosity). Both methods fail in reproducing this. They tend to stay in one region (the darker one in this particular case) and may lead to garbage growing. This is well represented with the synthesis maps where green is dominant for the example in row number five. Combining MSLG to [] used as post-processing significantly reduces this effect, although [] has its own limitations.

.. Influence of the parameters

In Influence of the patch size The synthesis results are very sensitive to the patch size, in particular for macro textures that have details at different scales. Figure . clearly shows that if the patch size is too small then the synthesis will fail. Using the multiscale approach strongly reduces the dependency of the method on the patch size. For all examples shown in this paper, the multiscale method using patches of 20 × 20 pixels was enough to guarantee correct synthesis results. For the one scale version the patch size used should have been much larger to produce convenient results, on the cost of reducing the variability of the Gaussian models and even creating verbatim copy effects. Influence of the overlap size For the nearest patches only the overlap areas are compared. This implies that the variance of the model estimated on that set of patches will be controlled only on the overlap region, thus allowing more variety in the remaining pixels of the patch. If that region is not big enough then the complementary region of the overlap will not be correctly modeled, since outliers can be considered in the set of patches. In

Influence of the number of neighbours



patches. More complex quilting techniques are then avoided. The multiscale algorithm fosters a better respect of the spread (but it's not a guarantee).

Input K = 1 K = 2 K = 3

. Conclusion

In this chapter, a local texture sampling method in the patch space using conditional Gaussian models was proposed. The motivation was to dispose of a patch stitching step by using a more robust local model for the texture. The Gaussian distribution of a patch was then conditioned to the values of its overlapping region. Two approaches were considered: CLG and RCLG models, and were compared to the local Gaussian model (LG). The results show that avoiding the stitching step is possible when dealing with periodic and pseudo-periodic textures. For more complex textures the conditional models are less performing and are fast limited by the size of the texture input sample. In general, the synthesis results using the conditional models were slightly smoother than the ones of the LG model.

The second contribution of this chapter was a multiscale texture synthesis algorithm using local Gaussian models (LG). In general, the experiments showed satisfying results for a wide variety of texture samples. Indeed, the patch-based approach conserves the local structures while sampling each patch from its Gaussian distribution, thus creating new patches that do not exist in the input sample. On the other hand the multiscale approach permits to synthesize the global arrangement of the salient structures of the input sample. The experiments also put in evidence that the use of Gaussian distributions, for some texture examples, have a tendency to slightly smooth the result  compared to the initial resolution of the sample. To overcome this effect the method was combined to Portilla and Simoncelli's algorithm [] as a first solution.

Several aspects of this method are still open to elucidation, such as the way patches are compared, the adaptation of the number of neighbours to the patch being modeled and the conservation of the texture's global statistics.



In this chapter we analyze the problem of stitching patches together. One of the major challenges is stitching patches that differ considerably along the edges. Stitching methods used in texture synthesis usually fail in these cases, resulting in evident contrast changes between the patches. We suggest a novel method based on optimal transport theory that manages to merge both patches producing a smooth transition between them. This approach works in two steps: first by finding a common edge using the midway equalization methods, followed by propagating this edge into both images producing a smooth transition. We thoroughly evaluated the proposed method first by trying it directly on selected patches, followed by integrating it in the texture synthesis method of Efros and Freeman. Based on the results the proposed method proved to be a possible alternative to perform patch stitching in nonparametric patch-based methods.

. Introduction

In this chapter we tackle the problem of stitching texture patches together. This problem arises repeatedly in the previous chapters where one of the steps of patch based texture synthesis methods is the fusion of overlapping patches. Given different patches of an input texture we want to stitch them together seamlessly. When the texture patches to stitch quite similar then existing techniques such as the quilting method in [], Poisson editing [], graph cuts algorithm [] among others give a satisfactory result. Nevertheless when the patches to stitch are not that similar on the blending edge then these approaches do not solve conveniently the stitching step. This can be observed in the illustration of Figure .. Here we address the question of deforming the content of the patches on a given curve in a way that both patches match on that curve. In other words we want to find a common signal such as the edges of the patches are transformed to match that signal. The problem is then composed of two steps. The first one consists in finding a common edge that is midway between the two edges. The second step consists in finding a way to progressively propagate this midway edge into both images. For the first step our method is strongly inspired in the midway equalization algorithm [].

.. Antecedents

The Poisson editing method The Poisson editing method [] is a method based on the manipulation of the image gradients for processing or combining images. It is used for seamless cloning, local illumination changes, seamless tiling etc.. In particular for seamless tiling the gra- dients of the images to tile are given as a vector field and the Poisson equation is used to find the image gradient field that is the closest to the given vector field. This technique allows to follow the variations of each region to tile while removing the visible edges between regions. The drawback we encountered using this technique for the purpose of tiling two texture patches is the fact that when propagating the seamless tile, the Poisson equation decays very fast and for some cases, even if the "real edge" disappears since the propagation is quickly attenuated a visual effect of edge still remains (see Figure .).

The quilting method

The quilting method presented in [] is a procedure to stitch a new patch in the sequentially built output texture by computing an optimal boundary cut between the patch and the synthesis area thanks to a linear programming optimization. This optimal boundary is the path on which the cumulative error between both patches is minimum. This technique presents very good results when the patches to tile are close to each other on their overlapping area. Nevertheless when this is not the case this optimal boundary is quite noticeable (see Figure .).

Adaptive stitching methods

Another work where texture merging is mentioned is in []. In this work Kwatra et al. propose to synthesize a new texture by copying irregularly shaped patches from the sample image into the output image. Once a candidate rectangular patch is selected an optimal portion of it is selected yielding an irregular shaped patch. The portion of the patch to copy is determined by using a graph cut algorithm, where the idea is to find out which patch each of the pixels in the overlapping area (of two patches A and B) come from. This amount to seek the minimum cost cut of the graph that separates the node A from the node B []. This method can work to merge patches. The method proposed in [] proposes to use a graph-cut formulation that uses α-expansions [] in order to take into account the previous seams and update them to deal with overlaps of more than two patches. They also propose an image merging technique using the same approach. Given a source image A and a target image B the idea is to constrain some pixels of image B where the image A will be inserted. Then the graph-cut algorithm finds the best seam that goes through the remaining unconstrained pixels of B setting the seam between both images. The results observed are satisfactory as long as both images have some common background. Otherwise, although the seam computed is optimal it will remain visible.

Our contribution

In this chapter we suggest an original method to merge two texture patches along a common curve (that we call the blending curve). The problem can be seen as finding a common one dimensional signal which is midway, up to a minimal deformation, between two one dimensional signals denoted s 1 and s 2 . The signals under construction are nothing but the restriction of the patches p 1 and p 2 to the blending curve. Afterwards the deformations sending s 1 and 2 to their midway signal are progressively propagated inside the corresponding patches. The signal registration step is inspired in the midway equalization method []. For this we considered that the two signals s 1 and s 2 were probability distributions on which we apply a geometric deformation derived from the midway equalization []. The common intermediary signal s is thus defined as the underlying histogram of the harmonic mean of the cumulative functions of s 1 and s 2 . We noticed that the behaviour of the registered signals using the harmonic mean of the cumulative functions not always gives as result the intuitive registration one could expect. The signals are well registered but the intermediary signal can sometimes be "unnatural" and this could complicate the propagation step. For the propagation step two alternatives are considered. The first one consists in linearly propagating the deformations. The second one considers convex combinations of two This chapter is organized as follows. Section . recalls the midway equalization algorithm []. In section . the midway blending technique is introduced. At first a description on a simple blending curve is given. It is followed by a generalization of the method. A description of the algorithm is provided as well as a discussion on the extension to a color version. As a final contribution of this section we describe how to integrate this method to a simple and efficient texture synthesis method. In section . experimental results are shown organized in two parts. The first one compares the midway blending method to [] and []. The second part compares synthesis results of [] to the same synthesis method using the midway blending. Conclusions are given in section ..



. Optimal transport -midway equalization algorithm

The theory of optimal transport [] has been widely used in computer vision as a metric between statistical features []. This theory takes into account the spatial location of the density modes and provides a transport plan between them. Depending on the probability densities that are being transported, the transport plans can be very irregular. In [] the authors propose a variational formulation to compute a regular transport map between two empirical densities. They show that the use of this regularized optimal transport improves visually the results for color image transfer.

The theory of optimal transport has been used in several applications, namely texture synthesis and texture mixing [, , , ], histogram equalization [, ], color normalization [] and color transfer [, , , ] among others. We follow by describing with more details the histogram equalization method [] which is strongly used in the remaining of the chapter.

Delon presented in [] an image equalization method. The comparison of two images is always better when both images have the same dynamic ranges and luminance. For this given two images I 1 and I 2 a contrast change is applied obtaining two new images with the same gray or colour intensity histogram. It is desired that the common histogram stands midway between the histograms of I 1 and I 2 staying as close as possible to both of them avoiding to follow one of them.

Let I : Ω → [0, 255] be an image and h the distribution of its values, i.e. the histogram of I.



The cumulative histogram of I is the increasing function H : [0, 255] → [0, 1] defined as

H(x) = x 0 h(t)dt. (.)
For a given function φ : [0, 255] → [0, 255] continuous and strictly increasing the cumulative histogram of the image φ(I) is H • φ -1 . Thus, applying the change of contrast φ 1 and φ 2 defined in (.) to the images I 1 and I 2 respectively transports the respective cumulative histograms H 1 and H 2 to a common one H. In the definition (.) H -1 is the inverse function of H. In all of the paper every instance of H -1 will refer to the inverse as a function.

φ i = H -1 • H i , i = 1, 2 (.)
In other words the principle of histogram equalization methods is to apply two contrast functions φ 1 and φ 2 as in (.) to the images I 1 and I 2 to bring I 1 and I 2 to have the same histogram H. There are several ways to determine this common histogram H. Among them the classical equalization method where H = id, the specification method where the common histogram is equal to one of the images cumulative histogram H = H 1 (or H = H 2 ), taking the mean of both histograms H = H 1 +H 2 2 . All the previously mentioned methods are in general unsatisfactory. In the first method it might be to drastic to flatten both histograms (if these are far from being flat), the quantization noise is enhanced with this change of contrast. The second method is suitable if the dynamic of both histograms is similar otherwise the one that is modified may try to follow inconsistently the dynamic of the specified image. The third mentioned method tends to create unexisting structures in the resulting images. In [] the author proposes to use the midway histogram H as the harmonic mean between the cumulative histograms H 1 and H 2 where

H -1 = H -1 1 + H -1 2 2 (.)
and then the contrast changes involved are defined as

φ 1 = id+H -1 2 •H 1 2 and φ 2 = id+H -1 1 •H 2 2 . (.)
In this way the equalized images have the same gray level distribution while preserving as much as possible their initial gray level dynamics. The midway histogram H is in fact the mid point of the geodesic in the Wassserstein distance [] between the distributions H 1 and H 2 .

. Midway blending

We extended the midway equalization method previously described to use it as a blending method for texture patches. For this we propose to deform the edges on which the patches are stitched in the spirit of []. We will consider the signals on these edges as histograms.

Let I 1 : Ω → [0, 255] and I 2 : Ω → [0, 255] be two images. Here to illustrate the blending method we consider that the images are blended along a horizontal edge as illustrated in Figure .. These edges will be referred to as the blending edges. Let s 1 : ω → [0, 255] be the one dimensional edge of I 1 and s 2 : ω → [0, 255] the one dimensional edge of I 2 with ω an interval [a, b]. We refer to h 1 : ω → [0, 1] and h 2 : ω → [0, 1] as the one dimensional histograms of s 1 and s 2 respectively. These histograms are nothing but the normalized one dimensional signals as defined in (.). The cumulative histograms H 1 and H 2 are then defined as in (.) for h 1 and h 2 respectively.

h i (x) = s i (x) b a s i (t)dt (.)
The deformations we aim at will expand and/or contract the two signals on the blending edges in order to remove the visible edge between I 1 and I 2 . We denote these two functions as φ 1 : ω → ω and φ 2 : ω → ω strictly increasing such as

H 1 • φ 1 = H 2 • φ 2 = H,
where the common cumulative histogram H is midway histogram defined in (.). The involved functions φ 1 and φ 2 are then defined as

φ 1 = I+H -1 2 •H 1 2 -1 and φ 2 = I+H -1 1 •H 2 2 -1
.

Let H1 and H2 be the deformed cumulative signals after applying φ 1 and φ 2 respectively

Hi (x) = H i (φ i (x)), i = 1, 2. (.)
where deriving (.) we obtain the deformed histograms

hi (x) = h i (φ i (x))φ ′ i (x), ∀x ∈ ω, i = 1, 2
V and the respective re-normalized deformed signals are

si (x) = mh i (x), ∀x ∈ ω, i = 1, 2
where m = m 1 +m 2 2 with m i = b a s i (t)dt for i = 1, 2. In this way s1 and s2 have the same mass.



Once the signals s 1 and s 2 are registered, i.e. the two images coincide on the blending edge, the deformations φ 1 and φ 2 have to be propagated. Otherwise we will create new discontinuities between the deformed signal s i , i = 1, 2 and its corresponding image I i , i = 1, 2. To avoid these inconsistencies we propagate the deformations φ 1 and φ 2 on a small part of the corresponding images I 1 and I 2 . This propagation is done progressively to attenuate gradually the deformation on a portion of the image that we call the propagation band (see Figure .). A first approach, that seems natural, is to use a linear propagation. That means that for each image I i , i = 1, 2 the midway deformation φ i , i = 1, 2 is gradually attenuated as defined in (.). Let Ω = [1, M ] × [1, N ] be the rectangular domain of the images I 1 and I 2 and to illustrate the blending edges let us suppose that s 1 (y) = I 1 (x 0 , y), y ∈ [1, N ] and s 2 (y) = I 1 (x 0 , y), y ∈ [1, N ] (see Figure .). We denote by W the width of the band on which we want to extend progressively the linear deformation. We denote by Ĩi , i = 1, 2 the images I i , i = 1, 2 after deformation which are defined as

Ĩi (x, y) = I i (x i , ỹi )J i (x, y), ∀(x, y) ∈ Ω, i = 1, 2 where J i (x, y) = f (x) + (1 -f (x))φ ′ i (y) and (x, ỹ) are the deformed coordinates applying Φ i (x, y) → (x i , ỹi ) = Φ i (x, y), xi = x, ỹi = f (x)y + (1 -f (x))φ i (y) (.)
The weight function f is defined as

f (x) = 1 if |x 0 -x| > W |x 0 -x| W if |x 0 -x| ≤ W .



Finally the images Ĩ1 and Ĩ2 are progressively re-normalized on the band on which the deformation was propagated to avoid any contrast change that may appear as an edge. The resulting images Ĩ1 and Ĩ2 are then

Ĩi (x, y) = Ĩi (x, y) m i (x) (f (x)m i (x) + (1 -f (x)) m) (.)
where m i (x) = b a I i (x, y)dy. Nevertheless, we notice that propagating linearly the midway deformation yields some visual artifacts that comes of the normalization. Indeed, the normalization step is important since we want to transport signals of the same mass. However, propagating linearly the midway deformation does not guarantee that the signals after deformation conserve their mass. This normalization artifact can be seen in the example of Figure .. A simple synthetic case is used to illustrate it where two binary images are to be blended. One can observe that the deformation achieved (in terms of shrinkage and dilatation) is coherent. Nevertheless, for the images Ĩ2 and Ĩ one can notice a shadowed area (marked with the red circles). This artifact in our opinion is the result of normalization step (see (.)). Thus a different normalization step should be considered if propagating linearly the midway deformations. The second propagation approach we considered for the midway deformations φ 1 and φ 2 consisted in a transition between the inverse of the cumulative histograms using a convex combination. We call it progressive midway. Let us suppose we applied the midway deformations to the signals s 1 and s 2 (blending edges of I 1 and I 2 ). Thus I 1 and I 2 match on their blending edge. Now for each image we aim at propagating the deformation on a small band (see A and H B 1 using a convex combination (see (.)).

H -1 δ (y) = (1 -δ)H -1 A + δH -1 B 1 , 0 ≤ δ ≤ 1 (.)
The histograms h m B 1 = b a I 1 (x 0 -W, y)dy.

s δ = h δ ((1 -δ)m A + δm B 1 ), 0 ≤ δ ≤ 1 (.)
This approach deals with the normalization problems evoked for the linear propagation case. However one can observe in the results of Section . that the deformations can sometimes be too abrupt. In the next section we analyze the behaviour of the midway algorithm as a blending approach. We will consider the progressive midway propagation approach.

.. Interpretation of the method

Let us consider two signals s 1 and s 2 having the same mass as illustrated in Figure . where s 2 is simply a translation of s 1 . We would expect the midway signal s mid between s 1 and s 2 to be their translations to their halfway point. Indeed, that is exactly the result of applying the midway deformation to both signals as can be observed in Figure .. As a second example let us consider another set of rectangular signals s 1 and s 2 , where the flat zones values' are greater than zero. One would expect to have the same behaviour as in the example of Figure .. However, the midway signal s is not a translation of s 1 and s 2 as one can observe in Figure .. The cumulative histograms show that transporting the mass of both inputs signals s 1 and s 2 creates an additional intermediary flat zone in s. This redistribution of the mass is not intuitive. In this simple case a solution would be to substract the minimum value to both signals s 1 and s 2 before their registration. This is possible here since s 1 and s 2 have the same minimum value. A third example is shown in We conclude from these experiments that applying the midway blending does not always yield an intuitive deformation to blend two images. In Section . we show some results on real texture images where one can observe that it works for several textures and fail for others. In general when blending texture images these "unnatural" behaviour are less evident.



.. Midway blending algorithm

For the algorithm description let us suppose that I 1 and I 2 are two discrete images defined on Ω = {1, . . . , M } × {1, . . . , N }. The blending edges are now defined as s 1 (j) = I 1 (x 0 , j), j ∈ {1, . . . , N } and s 2 (j) = I 2 (x 0 , j), j ∈ {1, . . . , N }. The respective histograms are then

h i (j) = s i (j) N k=1 s i (k) , j ∈ {1, . . . , N }, i = 1, 2 (.)
and the cumulative ones

H i (j) = j k=1 h i (k), j ∈ {1, . . . , N }, i = 1, 2. (.)
A summary of the method is given in Algorithm , Algorithm  and Algorithm .

Algorithm : Pseudo inverse

 Input: Function to invert H, x domain of H  Output: Inverse H -1 : k ← 2{re-sampling factor} : n ← length(H) {equal to length(x)} : for i = 1 to kn -1 do : z ← (i -1)/(kn -2) : j ← n l=1 1 (H[l]<=z)
:

Ih[i] ← 1 H[j+1]-H[j] ((z -H[j])x[j + 1] + (H[j + 1] -z)x[j]) : end for : for i = 1 to n do : H -1 [i] ← Ih[1 + (i -1)k] : end for : return H -1 Algorithm : Midway histogram  Input: cumulative histograms H 1 , H 2  Output: midway cumulative histogram H mid : n ← lengthH 1 {equal to length(H 2 )} : H -1 1 ← pseudo -inverse(H 1 , {1, . . . , n}) {see Algorithm } : H -1 2 ← pseudo -inverse(H 2 , {1, . . . , n}) {see Algorithm } : H -1 mid ← (H -1 1 + H -1 2 )/2 : H mid ← pseudo -inverse(H -1 mid , {0, 1/(n -
: m 1 ← n k=1 s 1 (k) : m 2 ← n k=1 s 2 (k) : h 1 ← s 1 /m 1 {d histogram of s 1 } : h 2 ← s 2 /m 2 {d histogram of s 2 } : H 1 ← zeros(n + 1) : H 2 ← zeros(n + 1) : for i = 2 to n + 1 do : H 1 [i] ← H 1 [i -1] + h 1 [i -1] : H 2 [i] ← H 2 [i -1] + h 2 [i -1] : end for : H mid ← midway(H 1 , H 2 ) {see Algorithm } : h mid ← zeros(n) : for i = 1 to n do : h mid [i] ← H mid [i + 1] -H mid [i] : end for : s mid ← h mid
: s 1 [j] ← I 1 [x 0 , j] {get blending edges from I 1 } : s 2 [j] ← I 2 [x 0 , j] {get blending edges from I 2 } : s mid ← midway -registration(s 1 , s 2 ) {see Algorithm } : s B 1 [j] ← I 1 [x 0 -W, j] : s B 2 [j] ← I 2 [x 0 + W, j] : Ĩ1 ← propagate(s mid , s B 1 ) {propagation step defined in Section .} : Ĩ2 ← propagate(s mid , s B 2 ) {propagation step defined in Section .} : Ĩ[i, j] ← Ĩ1 [i, j] if i ≤ x 0 Ĩ2 [i, j] if i > x 0 : return Ĩ .

. Generalization of the midway blending

Let I 1 , I 2 : Ω → [0, 255] be two images. We consider γ : ω → Ω a parametric curve defined as γ(t) = (x(t), y(t)), t ∈ ω and ω = [a, b]. Let s 1 : ω → [0, 255] and s 2 : ω → [0, 255] be the D blending edges defined as s i (t) = I i (γ(t)), i = 1, 2. Then the histograms h 1 and h 2 and the cumulative histograms H 1 and H 2 are defined as in (.) and (.). As explained in the previous  γ(t) For the orange curve the point (x(t), y(t)) moves (x(φ(t)), y(φ(t))). For the blue curve the point (x λ (t), y λ (t)) moves (x λ (φ(λ, t)), y λ (φ(λ, t))). For the gray curve the point does not move.

(x λ , y λ ) := W -λ λ γ(t) + λ W (γ(t) + W n(t)) γ(t) + W n(t) ψ(λ, t) γ(φ(t)) (x λ , ỹλ ) := (x λ (ψ(λ, t)) , x λ (ψ(λ, t))) γ(t) + W n(t)
section we aim at finding two functions φ 1 and φ 2 such as the new histograms h1 and h2 satisfy

hi (t) = h i (φ i (t))φ ′ i (t), i = 1, 2.
This implies then ∀t ∈ ω

Ĩi (γ(t)) = Ĩi (x(t), y(t)) = m m i I i (x(φ i (t)), y(φ i (t)))φ ′ i (t)
where m i = b a s i (t)dt, i = 1, 2 and m = m 1 + m 2 . We shall now consider all the curves γ α,λ , α ∈ {-1, 1}, λ ∈ [0, W ] that are defined on a band of width 2W around γ. These curves are taken in the normal direction n(t) = γ ′ (t) ⊥ γ ′ (t) and are defined as γ α,λ (t) = γ(t) + αλn(t).

We assume that the considered curves and the width W are such that ∀(α, λ, t 1 ) = (β, µ, t 2 ), γ α,λ (t 1 ) = γ β,µ (t 2 ) (see .). Once again we consider two propagation approaches: the linear propagation and the progressive midway propagation. The first one, the linear propagation (see Figure .), is defined as the deformation

ψ i : ω × [0, W ] → ω, i = 1, 2 ψ i (λ, t) = W -λ W φ i (t) + λ W t, λ ∈ [0, W ], t ∈ ω.
The resulting deformed images Ĩ1 and Ĩ2 are then

Ĩi (x α,λ (t), y α,λ (t)) = = I i (x α,λ (ψ i (λ, t)), y α,λ (ψ i (λ, t))) W -λ W φ ′ i (t) + λ W ,
for i = 1, 2, and after re-normalization

 Ĩi (x α,λ (t), y α,λ (t)) = Ĩi (x α,λ (t), y α,λ (t)) m i (α, λ) ( W -λ W ) m + λ W m i (α, λ) ,
where m i (α, λ) = ω I i (γ(u) + αλn(u))du.

Let us consider the particular case of a circular arc γ : ω → Ω. That is

γ(t) = (R cos θ(t), R sin θ(t)) = R(cos(s), sin(s)), t ∈ ω
where s = θ(t) = π 2 t, ω = [0, 1] and R is the radius of the circular curve. We consider the curve γ parametrized in the angles s. The normal n(s) to the curve γ(s) is defined as

n(s) = γ ′′ (s) γ ′′ (s) = - γ(s) R = -(cos(s), sin(s)) ,
where s denotes the arc length, i.e. γ ′ (s) = 1. The resulting images Ĩ1 and Ĩ2 are

Ĩi (R(cos(s), sin(s)) -αλ(cos(s), sin(s))) = Ĩi ((R -αλ)(cos(s), sin(s))) = K i (r)I i ((R -αλ)(cos(ψ i (s)), sin(ψ i (s)))) ψ ′ i (s),
where

K i (r) = ( W -λ W ) m + λ W m i (r) / (m i (r)) and m i (r) = π/2 0 I i (r(cos θ, sin θ))dθ.
We can verify that indeed for λ = 0 the deformation function is

φ i Ĩi (R(cos(s), sin(s))) = K i (R)I i (R(cos(φ i (s)), sin(φ i (s)))) φ ′ i (s)
and for λ = W, α ∈ {-1, 1} the deformation function is the identity

Ĩi ((R -αW )(cos(s), sin(s))) = K i (R -αW )I i ((R -αW )(cos(s), sin(s))) .
The second propagation method, the progressive midway, is as explained for the horizontal case: a convex combination of the inverse cumulative histograms of the blending edge and the second edge of the propagation area.

In practice we tested three curves: horizontal, vertical and circular arc. To simplify we transformed the vertical and circular case to a horizontal edge problem. A pseudo-code of the algorithm is provided in Algorithm . For this let us consider that I 1 and I 2 are two discrete images defined on Ω = {1, . . . , M } × {1, . . . , N }. The blending edges can be taken along a horizontal curve γ : {1, . . . , N } → Ω of length N defined as γ[j] = (x 0 , j) or along a vertical curve γ : {1, . . . , M } → Ω of length M defined as γ[j] = (j, x 0 ) or along a circular curve γ : {1, . . . , L} → Ω of length L defined as γ[j] = (Mx 0 cos θ[j], Nx 0 sin θ[j]) where θ[j] = j-1 L-1 π 2 , j ∈ {1, . . . , L} and x 0 is the radius of the curve. The blending edges are then defined as I 1 ← T (I 1 ) {transform to a horizontal edge case}

s i [j] = I i (γ[j]) i = 1, 2.
:
I 2 ← T (I 2 ) {transform to a horizontal edge case} :

x 0 ← t(x 0 ) {transform to the reference position of a horizontal edge case}

: end if : Ĩ ← midway -blending(I 1 , I 2 , W, x 0 ) {see Algorithm } : if edge-type = horizontal then : Ĩ ← T -1 ( Ĩ) {transform to the original edge case} : end if : return Ĩ ..

 Colour midway blending

The first approach we considered for colour images was doing the transport on the grayscale image and then applying the same deformation to each channel. Thus no color artifact appears. Nevertheless for each colour channel the signals of both images do not match perfectly and an edge is still evident on the resulting colour image.

As a second alternative we considered treating each channel independently. This amounts to apply a midway blending deformation for each channel separately. In this case the edges do match perfectly but false colors may appear, each colour channel being processed independently from the other ones. Indeed, the colour components of the RGB color representation are strongly correlated. Thus, deforming each channel independently will not maintain this correlation. The solution is to work with a different color space where the components are less correlated.

Instead of working with the RGB color spae we propose to use a DCT colour representation which is similar to the HSV colour space but easier to implement. Let us suppose a colour image I : Ω → [0, 255] 3 defined in the RGB space. Than the DCT representation of I is obtained from the RGB components applying the orthogonal linear map

T =   1/ √ 3 1/ √ 3 1/ √ 3 1/ √ 2 0 -1/ √ 2 1/ √ 6 -2/ √ 6 1/ √ 6   .
Notice that the DCT representation of I yields negative values in its second and third components (the first component is the average of the R, G and B channels up to a constant multiplier). Since the signals are considered as histograms and their cumulative functions are used, negative values are not desired. To avoid this after applying the orthogonal matrix T we substract its minimum from each channel, so that it has no negative values. Another approach one could consider is the PCA colour space where the colour representation is not fixed and depends on the image. This approach yields a representation where the colour components are decorrelated. However when using a PCA colour space each image has a different diagonalization basis and the question that arises from this is what common basis should one use to apply the transport between signals? We observed that for natural images the DCT colour  representation yields components that are in average almost decorrelated. This representation is very close to the PCA colour space with the advantage that the same basis is used for all images. Nevertheless to our application the patches we aim at stitching come from the same texture image. Thus applying the PCA decomposition to both images yields an orthogonal transformation very similar to the one that we would obtain if the PCA decomposition was applied to each image separately. We compare the three approaches (RGB, DCT and PCA) in

.. Texture synthesis

This work is originally motivated by patch based texture synthesis methods. Up to now several stitching techniques were used, among them the quilting method [], graph cuts [], Poisson blending [], etc. However, these techniques do not yield an image deformation. They find the "best way" of putting together the different pieces. Methods [] and [] tend to find the best possible path in a given region without changing the content of the image. For instance in the example of Figure . the stripes are not placed one in front of the other. The visual impression of an existing edge between both images remains. A different approach would be to use the method to seamlessly edit image regions [] to blend the patches. Still, one can notice that the edge is attenuated but remains evident and the stripes do not match. In this chapter we aimed at "matching these stripes". To that end we allow to stretch and/or dilate the images on their common edge as illustrated in Figure . where one can observe how the stripes "moved" and match.

In this section we detail the synthesis algorithm where the midway blending method is used. We shall compare three stitching techniques: the Poisson editing methods [], the quilting method in [] and the midway blending. In continuation, we shall compare two texture synthesis methods: one using the quilting technique [] and another using the midway blending.

The image quilting method [] computes a new texture by rearranging seamlessly patches of the input texture. The procedure to stitch a new patch in the sequentially built output texture is achieved by computing an optimal boundary cut between the patch and the synthesis area thanks to a linear programming optimization. This method is a sequential patch-based algorithm. I and Ĩ denote the input and the output texture respectively. The output image Ĩ is constructed patch by patch in a raster scan order. The goal of each iteration is to fill a patch of Ĩ that is only partially defined on a region called overlap region. Note that there are three possible overlap regions: vertical overlap for the first row, horizontal overlap for the first column, and L-shaped  overlap everywhere else. To complete the patch P A a patch P B of I that matches P A on the overlap region is randomly selected. An optimal boundary cut between P A and P B is then computed within the overlap region. This optimal boundary cut is used to construct the new patch P by blending P A and P B along the cut. We now want to compare the results of this texture synthesis method when replacing the stitching step by the midway blending. For details on the quilting method please refer to Annex A... The synthesis algorithm using the midway blending is detailed in Algorithm . Select a compatible patch P B ∈ I :

Construct P blending P A and P B using the midway blending algorithm (see Algorithm ). Depending on the overlap type the (blending curve) can be a horizontal, vertical or circular curve. (see Algorithm )

:

Replace P A with P within Ĩ :

return Ĩ : end for

. Experiments

.. Comparison of blending techniques

In this section we compare the blending methods that can be seen in Figures . to .. Given two patches p 1 and p 2 we applied three different stitching techniques: the quilting method [],

Poisson blending [] and the midway blending. For each texture example we blended the two patches along three different curves: horizontal, vertical and circular. For the quilting method the circular curve corresponds to the l-shape overlap. To compare the circular curve and the Lshape overlap we consider the same overlap size and the circular curve is taken such as its radius is contained in the l-shape overlap. For the midway blending results the propagation step used is the linear propagation approach. We noticed that the linear propagation yields better results for natural images and the normalization artifact is not evident as in Figure ..

For the horizontal and vertical blending edges one can observe in general that the results are comparable to the other methods and sometimes better when we can observe the quilting path or where the Poisson blending "smoothes" the edge but this one is still evident. However some results using the midway blending may present deformations that are too abrupt ( These failures are colour artifacts or highly abrupt deformations that are too evident. On the other hand, in the mentioned figures we attempted blending of strongly dissimilar patches, so that fusing them is close to a mission impossible. For  the circular edges however the results are not as convincing as for horizontal or vertical edges. This will depend on the curve γ used and the discretization of it. Nevertheless one can observe that when used in patch-based texture synthesis since the patches are quite similar on the overlap zones (it is intended to be this way) then using a circular curve is not problematic. One can also notice fewer or less color artifacts for the same reason.

.. Comparison of texture synthesis results

In this section we evaluate the results of synthesizing textures with patch based approaches using as stitching step the midway blending algorithm. We compare the results to [] in Figures . to .. The synthesis algorithm is the same in both cases, only the stitching step is modified. To compare the results the synthesized images are both initialized with the same seed patch and the parameters are exactly the same (patch size n and overlap size o).

In general the results of both methods yield satisfying and comparable visual results. Thus the midway blending approach could be a satisfying alternative to stitch patches together. This is particularly so for textures having a dominant structure direction, as for example the wood texture in Figure ..

An interesting aspect of the stitching method would be to relax the constraint on the resemblance of the patches on the overlap area. This could be interesting since as pointed out in Annex A. when relaxing the error tolerance on the overlap resemblance of the patches more innovation is introduced in the result. Increasing this tolerance error value permits to increase the number of candidate patches used at an iteration of the synthesis method and so the chosen patch may vary more. As a consequence the patches P old and P in might be more distant to each other on the overlap area and so by just using the quilting method to stitch them the boundary might be noticeable since the values of thoses pixel are not alterated. Nevertheless with the midway technique a seamless fusion remains possible.

. Conclusion

In this chapter we have presented a method for stitching texture patches inspired in the midway equalization algorithm []. This method consists in equalizing two images by transporting their gray levels cumulative histograms to a common one, the midway histogram. The cumulative midway histogram is the harmonic mean of the cumulative histograms of both images. The midway blending algorithm proceeds in the same way. Thus we consider the two signals we aim at registering as two distributions. The signals s 1 and s 2 on the blending edge are handled as histograms. The registered signals coincide with the midway signal defined as the histogram whose cumulative function is the harmonic mean of the cumulative functions of s 1 and s 2 . The midway blending is composed of two steps. The first one consists in registering the signals obtained by restricting both patches to their common boundary. The second one consists in propagating the two midway deformations toward the interior of both patches. So far we have observed that the registration step achieves satisfying results whereas the propagation step is sometimes more challenging. We explored two approaches for this. The first one, the linear propagation, sometimes yields normalization artifacts. The second approach, the progressive midway, has no normalization issues. Nevertheless applying the second approach for natural images is too abrupt and the deformations are sometimes too evident. However for the linear propagation the results on natural images are satisfying and the normalization artifacts are almost invisible, unlike what we observed for the  synthetic image examples. The examples in section .. were designed to illustrate that using the harmonic mean of the cumulative histograms does not always yield an intuitive solution.

In the first part of Section . we showed the result of blending two patches of the same texture by different stitching techniques. We observed our results are generally comparable or better. The quilting method [] does not modify the content of the image. Instead, it looks for the edge where the patches best match. The Poisson editing method [] yields a solution which mitigates the edges. Still, in challenging cases, both methods can leave behind conspicuous stitching edges. The midway blending avoids such edges, but it may sometimes require a patch deformation that is too abrupt. Besides color artifacts can be observed which are not present with the other techniques. In the second experimental section we tried the proposed technique in a patch based texture synthesis algorithm [] where we replaced the quilting step by the midway blending. We observed that the results obtained with the midway blending are comparable to those obtained with the quilting technique. Since for such algorithms the patches to blend are similar on the overlapping area the artifacts caused by the attempts at blending too dissimilar patches no longer occur. We conclude that our technique could be used as a stitching step for this kind of texture synthesis algorithms.

The method as we have seen has some limitations, such as the normalization issue for the linear propagation. Other techniques could also be considered with the aim at stretching and/or expanding the images on the overlap area. We showed that the proposed method can be use as a stitching step in patch based texture synthesis methods and has room for improvement as for more complex textures to blend. It is clear at this point that the deformation technique we introduced should be explored further, and perhaps as generalized, by using other deformation and registration tools such as optical flow. In this chapter we aim at analyzing the limits of exemplar based texture methods when using large natural texture images, a new problem which we might call Big Exemplar Based texture synthesis. Unlike exemplar-based methods that use small texture samples the goal is to synthesize large textures. We shall explore this problem with the tools at reach, namely the best methods from the extensive literature, and will not hesitate to combine them. To anticipate our conclusions, none of the existing methods manage by itself to satisfactorily emulate "big textures". Yet, we'll find that a clever combination of these methods yields interesting-if not perfect-results. Interestingly, the existing methods provide very good results on small crops of these images. But we discovered with some dismay and shall illustrate on many examples that "big textures" are very different from academic texture crops. They show drastic local textural variations, even if their overall effect still is the one of a single texture. This study is therefore a good tested to put on trial all existing exemplar-based texture synthesis methods.

















. Introduction

This chapter shows an extensive study on the limits of exemplar-based methods when used to synthesize large natural textures. The goal of these methods is to create from a given sample decor (for example a wood texture as in Figure .) new, different samples that are perceptually equivalent to the input example. Current exemplar-based approaches use small input samples. When the texture example is not stationary, non-parametric methods often end up in a garbage growing, meaning that only a local piece of the texture is being reproduced to the detriment of the rest. The parametric methods present a different problem in presence of a non-stationary texture: when the input sample contains more than one texture global statistics are learnt on the input. It follows that the result is a homogenization of the input sample. The texture images we try to synthesize in this chapter are quite complex: large and non stationary.

We began by trying some of the exemplar based texture synthesis approaches (Efros and Freeman [], Portilla and Simoncelli [], Galerne et al. [] and the multiscale locally Gaussian method presented in Chapter ) on the following type of images: wood, stone and metal (Figures ., ., ., ., . and .). Among these textures the most challenging were the wood samples. We observed that the synthesis results for the wood images were not satisfying. However for local crops of these textures where no strong global organization were evident the results obtained are impressive. One concludes obviously that, in spite of their perceptual unity, large textured images are more than a texture in the usual academic sense. In Figures .  are shown where one can clearly see that within a same texture very different textures are present. Wood planks for example are non-stationary images, containing several different textures and with very recognizable and unique structures such as ring and knots. For other complex surfaces such as the stones and metals the exemplar based techniques seem to give more satisfactory results. Again these textures are not as the academic textures considered for the experiments done to evaluate their performance. These images are more complex as can be seen in Figures . All the experiments done are showed in the following sections, including failures that permit to discard certain techniques for certain texture. We follow by showing results of a combination of the considered texture synthesis methods, there are two of them:

• The Efros and Freeman method [] combined to the locally Gaussian method

• The Portilla and Simoncelli method [] combined to the locally Gaussian method Combining the methods yields better visual results in particular for the combination PS + LG. Still this combination will not make it for big textures. We then show other emulation attempts, like for example the application of geometric deformations to the inputs aiming at changing the aspect of the salient structures. Randomly switching sub-parts of the input to partition the salient structures and reorganize them in a way that in the outputs the repeated patterns are less evident. As a last experiment we tried to vary the low frequency information while leaving unchanged the high frequencies with the same purpose. In all three emulation attempts we observed interesting behaviours for some of the textures.

This chapter is organized as follows. Section . describes the type of inputs we worked on. In section . a description of some of the state of the art exemplar-based texture synthesis methods is provided listing the pros and cons of each of them. In section . a first set of experiments are showed on small crops of the big textures. All of the exemplar-based methods of section . are tested on these images. They are tested separately and then some of them are combined. The results on these images are very satisfying. The same experiments are performed on the big texture themselves to try to emulate the structures at low resolution. The results are not very satisfying since the typical drawbacks of exemplar based methods arise (verbatim copies and garbage growing). This brings us to test other alternatives presented in Sections ., . and .. We try the application of one and two dimensional anamorphosis to the inputs, a random combination of slices of the input texture and the application of the [] method to randomize the low frequencies of the input.

. Large natural textures

Wood, stone and metal textures In this section we present images of typical "big textures". These images are of very high resolution(more than 1GB). For this reason the experiments were sometimes done on zoomed out versions of these inputs. We considered three types of classic natural textures: wood, stone and metals shown in Figures ., ., ., ., .. As a general remark on these inputs one can notice that they hardly represent "real textures" in the academic sense. Indeed we observe on different crops of the inputs, of size 500 × 500 pixels, that within one of these large natural textures several different textures are present (Figures .



provided show four completely different textures. Still worse, these crops themselves are not stationary. This obviously illustrates how hard it will be to synthesize a new different textures from such an example.

Among the inputs that we considered the more challenging ones are the wood textures. Strong salient configurations are evident in these images that makes difficult to generate new variants of these structures instead of mere "copies" of them. The stone and metal decors are in this sense easier. They do not present strong salient patterns. Nevertheless they still are highly non stationary.

. Exemplar based texture synthesis methods

Exemplar-based texture synthesis is defined as the process of generating from an input texture sample a perceptually equivalent larger one. Exemplar-based texture synthesis algorithms are commonly divided into two categories, the statistics-based methods and the non-parametric patchbased methods. The first category models a given texture sample by estimating statistical parameters that characterize the underlying stochastic process of the input. Although these methods can faithfully reproduce some of the global statistics of the sample and synthesize micro and pseudoperiodic textures, they generally do not yield high quality visual results for more structured ones, in particular when the sample is small and contains large objects. The second category rearranges local neighbourhoods of the input sample in a consistent way. These methods provide efficient algorithms able to reproduce highly structured textures. Even though they yield visual satisfactory results, they often turn into practising verbatim copies of large parts of the input sample and growing garbage.



Within these two groups we chose some of the state of the art methods to synthesize the given decors in section .. These methods are:

• Portilla & Simoncelli [], denoted PS • Galerne et al. [], denoted RPN • Efros & Freeman [], denoted EF • Multiscale locally Gaussian (Chapter ), denoted MSLG
We follow by giving a short explanation of these methods.

.. The Portilla and Simoncelli algorithm

The Portilla ans Simoncelli method [] is performed in two steps: analysis and synthesis. First, the texture sample is characterized by a set of statistics estimated on the coefficients of its steerable wavelet decomposition []. The steerable pyramid is an overcomplete linear multi-scale and multi-orientation image decomposition. The statistics extracted on the coefficients of this decomposition are their autocorrelations, cross-correlations (inner and intra scales) as well as the statistical moments of order one, two three and four of the input sample's values. Second, these statistics are enforced on a white noise image obtaining like this a texture satisfying these constraints. To this end the authors in [] suggest to project the synthesis image into the subspace of constraints iteratively until a stabilization of the image observed. The final image may not have exactly the same statistics as the input sample but reached a local minima where it stabilized. The synthesis results in a pre-attentive examination are in general indistinct from the original samples. Nevertheless to an attentive examination the synthesis results on structured textures often present blur and phantom effects. The structures are missed and the method tends to homogenize the output texture. We follow by showing in



.. Micro-texture synthesis by phase randomization

The random phase method [] synthesizes a texture by randomizing the phase of the Fourier coefficients of the input sample. The results are very satisfying for textures that are characterized by their Fourier modulus. This method will not grow garbage or verbatim copies of the input sample. It is however limited to micro textures. It will fail synthesizing structured textures. We follow by showing in Clearly, the knowledge of the modulus of the Fourier coefficients of these textures is not sufficient to recover the strong contrast of the inputs. For a detailed description of the method please refer to [].

.. The Efros and Leung algorithm

The Efros and Leung algorithm relies on Claude Shannons's Markov random field model for the English language []. He proposed to use a Markov chain to model the underlying stochastic process of a sequence of letters in a piece of English text. The same idea was used by Efros and Leung in [] to synthesize a texture image by considering that a pixel value depends on the values of its neighbouring pixels. The method is illustrated in Figure . and works as follows. For a given input texture a new image is synthesized pixelwise. For a pixel p being synthesized the algorithm finds all the neighbourhoods in the input image that are similar to the neighbourhood of p. Then one of this neighbourhoods is randomly chosen and the central pixel value is affected to the pixel being synthesized. The neighbourhood of p is a square patch centered in p and only the known pixels of this patch are considered when comparing to the neighbourhoods of the input. In general the visual results are very impressive most of all for structured textures. Nevertheless this algorithm suffers from two important drawbacks: verbatim of the input and garbage growing (the algorithm starts reproducing iteratively one part of the example and neglects the rest). We invite the reader to refer [] for a detailed description of the method.



.. The Efros and Freeman algorithm

Efros and Freeman's method [] is an extension of Efros and Leung's method. It is based on the same principle where the pixel values are conditioned to their neighbourhood values. Efros and Freman propose to generate a new image patch by patch in a raster scan order. At each iteration a patch that is only partially defined on a region called overlap region is completed. This is the patch under construction. To do so a patch of the input image that matches the patch under construction on the overlap region is randomly selected (patch selection step). An optimal boundary cut between the chosen patch and the one under construction is then computed within the overlap region (stitching step). This optimal boundary cut is used to construct the new patch by blending the chosen patch and the patch under construction along the cut. The results obtained are very impressive, particularly for highly structured textures. The patch size being larger, the risk of garbage growing is reduced comparing to the Efros-Leung algorithm. Nevertheless the risk of verbatim copies remains ans is even amplified. Moreover, the respect of the global statistics of the input is not guaranteed and this is quite visible when the input texture is not stationary (for example if there is a change of illumination across the image). In This follows the fact that the input texture is not homogeneously illuminated. However for the second example where this illumination effect is not evident the synthesis is excellent. The last two



.. The multiscale locally Gaussian algorithm

The multiscale locally Gaussian method presented in Chapter  is a parametric patch based method. The output image is constructed patch by patch, in a raster scan order. In contrast with the Efros and Freeman work [] this method involves a multiscale local Gaussian model in the patch space. Every patch used to construct the output texture is sampled form a multivariate Gaussian distribution (a different one for each patch). For a given patch p under construction the set of m patches from the input texture that are similar (to the one under construction) to p on their overlap area is used to estimate the parameters of the multivariate Gaussian distribution used to sample a new patch. This is done for each patch in the output image. Figure . illustrates an iteration of the method. The multiscale approach permits to use "different patch sizes" within one synthesis and hence capture the details at different scales. The results are visually satisfying. The quality of the results is higher than those of statistics based methods. Nevertheless, even though the patches generated do not exist in the input texture the risk of visual verbatim copy remains.

. Applying the classical methods

.. Experiment : exemplar-based methods

We started by applying the exemplar-based methods described in section . to small crops of the provided inputs. In In general when applying RPN and PS to the small crops the results obtained are too blurry. In the case of EF and MSLG, for some inputs, satisfying results are observed. This is true especially when the input is stationary. Otherwise garbage growing and verbatim copies are often a risk. One can also observe that for some examples the result of the MSLG method might be slightly blurry. In Figure ?? larger crops of the inputs are considered. These crops are taken on images of lower resolution regarding the original inputs. For these examples it is evident that applying the PS or RPN method will not work at all. On the other hand for the EF and LG algorithms the verbatim copy and garbage growing is more evident and likely to happen, particularly when the inputs have strong isolated patterns, as it is the case with wood textures. We follow by showing a combination of these methods on the same examples.



.. Experiment : combining exemplar-based methods

In this section we present some synthesis results when using a combination of different texture synthesis methods. We apply them to the two sized crops of the different input images. The combinations considered are the following:

• The locally Gaussian method combined to the Efros and Freeman method denoted as LG+EF • The locally Gaussian method combined to the Portilla and Simoncelli method denoted as

LG+PS

LG + EF The combination of the locally Gaussian method with the Efros and Freeman method consists of two steps. The first step synthesizes the given input I with the locally Gaussian method generating a new texture image that we denote I lg . The second step consists in applying the Efros and Freeman algorithm to the given input sample initializing the output image that we denote I ef with the image I lg . The method is basically the same as the one described in section .. The only step of the algorithm that is modified is the patch selection step. In the method described in [] at each iteration the added patch was chosen among those (in the input sample) whose overlap region was similar to the one of the patch under construction. When combining the methods instead of only comparing the overlap areas we compare the entire patches. Initializing the output with a first synthesis I lg enables the method to use the whole patch under construction to find a candidate in the input sample I. The candidate patch taken from I is then quilted in I ef in the corresponding position with the same stitching step as in []. This combination allows to recover the lost resolution of the LG synthesis. However it is not capable of masking the garbage growing effects as effectively LG+PS combination does.



LG + PS The combination of the locally Gaussian method with the Portilla and Simoncelli method consists of two steps. A first step where given the input image I a new texture I lg is generated using LG. The second step uses PS where the initialization "noise image" is replaced by I lg . As explained in section . the statistics to impose are learnt on the input I. What follows is a synthesis step where the output image is projected on the subspaces of constraints. There exist several local solutions to this projection step. When initializing PS with the result of LG the initialization image is generally quite close to the images living in the sub-space of the whole set of constraints. Thus the result obtained is improved compared to PS images starting from a random noise image. Naturally fixing the initialization of the PS algorithm removes the randomness of the generated texture. But this is not the case since the initialization is itself random as it is generated from another random process. We observed that combining EF to LG sometimes slightly improve the granularity of the result but that if the result of the LG method diverged then applying EF would not mask these problems. However, the combination of PS with LG masks efficiently these drawbacks ans improves significantly the quality of the synthesized image.

In Figures . and . one can observe that so far the combination LG+PS permits to emulate all local textures with no definite global structure. That is for local synthesis combining these two methods (LG and PS) is a good solution to synthesize locally the decors provided. These small images are easy textures. We now move one step farther and attempt a synthesis of larger parts of the original inputs (see Figure .). As expected the results are less satisfying than for local synthesis.

In summary, LG+PS is a sophisticated method, more than any existing, and it works satisfactorily for local textures with no strong global organization. Thus it makes correct and even amazingly good imitations of small pieces of textures, yet remains unsatisfactory for large structured pieces of wood. Thus another method must be considered for such images which in fact are more than a  texture, as they have this non-stationary global organization with rings and knots, among other special, very recognizable structures. In short, wood at a large scale is NOT a texture in the usual sense and therefore imitations of it will either fail or tend to repeat the input. We therefore tried to expand the method to capture a structure at low scale and scale up with details that we know how to do deal with it. As we shall see, this will fail because it will fatally create repetitions.

.. Experiment : multiscale on the low resolution inputs

A last experiment was performed with the locally Gaussian algorithm using the multiscale approach (MSLG). The goal was to create a low resolution landscape on which details are progressively added. The multiscale approach is a classic solution to contemplate details at different scales within an image. This approach was applied on a low resolution version of the whole sample images.

The existing methods are able to perform a correct synthesis on some details of the provided examples at several scales. However due to the non stationarity of the provided textures, adaptations of the methods should be considered. We observed that even using a multiscale approach to create a low resolution landscape is not enough for the wood texture where large and special global organizations are visible and difficult to synthesize without creating a verbatim copy effect or garbage growing due to the strong non stationarity of this kind of images. Hence we decided to see if some perturbations of the input might be enough to hide the resemblance applying some geometric deformations on the inputs. This is presented in the following section.



Figure 5.19 -Overview of the multivariate Gaussian distribution of a patch. For the patch q (pink patch in the output image) being synthesized the algorithm finds the m nearest neighbours in the input image that matches q. Then, on this set of m patches a multivariate Gaussian distribution is estimated from which a patch p is sampled and quilted in the output at the corresponding position.

. Masking repetitions with anamorphosis

As we have seen in the previous three experiments, applying and combining different exemplar based texture synthesis methods works for local textures. For more complex images this is not the case. The images are highly non stationary and this makes it very difficult. Even though we tried synthesizing the inputs at lower resolution aiming at recovering a synthesized sample containing the global arrangements to continue by adding progressively the details of the different scales, the results were not satisfactory. The textures still retain a strong global organization at lower scales. If the low resolution synthesized image contains some repetitions, for example, then it is not easy to recover when scaling up in resolution. Thus we decided to change the strategy and aim at hiding these undue repetitions. As a first alternative we considered applying different anamorphoses to the inputs in order to deform the global organizations of the sample in a way that it is unrecognizable from one synthesis to another. Yet, the deformation cannot be too brutal. For example in a wood image, the deformation of a knot can be strong but should maintain it approximately convex. Unfortunately, we observed that after deformation, these global configurations are still recognizable from one image to the other, in particular for wood textures. We tried several different deformation functions that we separate in one dimensional anamorphoses and two dimensional anamorphoses. Basically, the one dimensional anamorphosis will deform the input in one direction and the two dimensional anamorphosis in two directions. For the one dimensional case we chose three simple functions to analyze their behaviours. For the two dimensional case we chose a geometric deformation performed as a simulation of turbulence [].

Applying several different one dimensional and two dimensional anamorphosis to the wood





. Random mosaics

Another approach we considered to avoid evident repetitions between synthesis and within one same synthesized image is the random mosaic method. This experiment consists in slicing the input image vertically (or horizontally), randomly permuting these slices and after blending together the pieces using a stitching step to recompose a new texture image. The stitching technique can for example use the approach of image editing in []. The advantage of this approach is that no "unnatural" repetition of the global configurations will be observed since each part is used once. We also hope that the slicing succeed to separate the salient structures but not up to the point where the salient and less salient structures would be completely lost. This approach restricts the new sample to have the same size as the input and the number of possible permutations is limited.

A possibility to overcome these limitations was to consider more than one input image in order to generate one new texture. There is no reason to limit ourselves to the use of a single input sample as it is done up to now in the exemplar-based methodologies.

The following experiences show three variants of this approach. The first one called the one dimensional random mosaic consists in slicing the image in one direction direction (horizontal or vertical). The second variant, the two dimensional random mosaic, consists in slicing the input in both directions, thus generating patches that are put together after being randomly permuted. This variant is suitable for textures that are isotropic like the stones and metal images that we presented. For the wood textures this variant is not of interest since there is an evident dominant direction and partitioning the image in both directions would completely alter the nature of the input image. The last variant, the multi-random mosaic, consists in using several inputs and then applying the one dimensional or two dimensional random mosaic, randomly permuting all the pieces and constructing a new images by picking some of the pieces.

Applying the random mosaic algorithm yields very good results for the case of the stone and metal textures. The generated images are all very different from one synthesis to another. Controlling the size of the slicing is important in order not to destroy the characteristic structures of these images. However for the wood texture this approach is not suitable. All global organization are lost and some of the salient structures are recognizable from one synthesis to another. For the case where the random mosaic approach works one can as well observe that the generated textures look different at all scales seeing the textures from far away still yields visual results that are different.

. Synthesizing low frequencies

One of the difficulties of the input images given in section . is to correctly synthesize their low frequencies, i.e. their global organizations. We have seen with the previous experiments that using a multiscale approach is not satisfactory, since repetitions are still visible at lower scales hence propagated to the higher scales. Ideally we aim at generating a low resolution image having the "natural" amount of global configurations on which we gradually add the correct granularity of the different scales.

For this purpose we suggest to separate the input image I into two components: the low and high frequency images. We denote by G θ a Gaussian kernel and define as the low frequency image I LF = I * G θ and as the high frequency component I HF = I -I * G θ . We follow by applying the RPN method [] to the low frequency component I LF . In other terms, we add a random phase σ to the Fourier transform ÎRPN ← ÎLF e iσ where σ is a white noise image niformly distributed over [-π, π]. This new low frequency image is denoted as ĨRPN . To recompose the new texture we simply cumulate the low and high frequencies components Ĩ = ĨLF + I HF .



. Mixing texture images

The last experiment that we considered consisted in mixing two texture images. Ideally one of the textures denoted as I 1 should be very sober and the second texture I 2 would have no restriction. Let us consider the case of wood textures. This could allow us to dispose of several sober images of a same essence of wood on which we will add different particular structures taken from a second image. This approach is clearly limited with respect to the number of images one can generate. This approach is considered as a tool which allows us to add some variations on the generated textures. To achieve the texture mixing we used the Poisson editing work described in [].

The goal of the Poisson editing technique [] is to seamlessly edit and clone selected regions of an image. The mathematical tool at the heart of the approach is the Poisson partial differential equation with Dirichlet boundary conditions. Given two images I 1 and I 2 , the goal is to replace a region of I 1 with one of I 2 . The domain to replace in I 1 is denoted by Ω and δΩ its boundary. Then the solution that we denote as I minimizes the following equation

∇I(x) = div( I 2 ) if x ∈ Ω I(x) = I 1 (x) if x ∈ δΩ.
Another application of the Poisson editing method [] is to add or remove details to an image mixing their gradients retaining the stronger variation at each pixel of the image. This is what we use to mix two different texture images obtained by applying (.). The results obtained are interesting. However the synthesized images are still quite similar to the input ones.

∇I(x) = ∇I 1 if ∇I 1 ≥ ∇I 2 ∇I 1 if otherwise. (.) In Figures .

. Conclusion

In this chapter we have presented an extensive analysis on how to imitate large natural textures and explored the limitations of this endeavour. The texture images considered in this chapter were high resolution photographies of different materials such as wood, metals and stones. The study shows that, in spite of the remarkable progress that have been done in the area of exemplar-based texture synthesis methods, they are posed for a very particular problem that is using stationary and small samples of image textures. The existing methods are still incapable of dealing with such complex textures containing strong global configurations and being non-stationary.

We have tested several methods and combination of methods. In general the results achieved on small crops of the given input images are satisfying. In particular, the combination of the locally Gaussian method with the Portilla and Simoncelli method yields satisfying results for small images. For bigger images presenting strong global organizations all methods fail. They cannot work on highly structured and non-stationary images. Insisting on using such exemplar-based methods leads to repetitions of the conspicuous structures.

We tried other alternatives such as applying random permutations on different parts of the input image or applying one dimensional and two dimensional geometric deformations to mask the conspicuous structures but none of these successfully hide the repetitions to the human eye. This conclusion is final in the case of the wood textures who have very strong salient configurations. The one dimensional and two dimensional random mosaic method returns satisfactory results on some samples.



As a more general conclusion when trying to synthesize samples that are not stationary and very large the exemplar-based methods fail. Patch based methods yield too many repetitions of the salient structures. The risk of growing garbage exists due to the strong non-stationarity. Methods like the Portilla and Simoncelli are limited by the large size of the input. Combining the Portilla and Simoncelli algorithm with the locally Gaussian method works, up to a given resolution. Afterwards the methods tend to homogenize the input; i.e. will start "mixing" the different structures present destroying them. We deduced that using a single wood texture sample, it is almost impossible to generate several different new images from this example. It is very difficult to randomly create knots with just one example without repeating it. This is a big limitation for exemplar based methods. The exemplar based methods usually work because the examples images present a lot of auto similarities that can be used on the one hand to infer a stochastic process defining the texture or on the other hand to estimate the models of a patch-based method.



A Efros and Freeman Image Quilting

Algorithm for Texture Synthesis

Exemplar-based texture synthesis is defined as the process of generating, from an input texture sample, new texture images that are perceptually equivalent to the input. Efros and Freeman's method is a non-parametric patch-based method which computes an output texture image by quilting together patches taken from the input sample. The main innovation of their work relies in the stitching technique which significantly reduces the transition effect between patches. In this paper, we propose a detailed analysis and implementation of their work. We provide a complete mathematical description of the linear programing problem used for the quilting step as well as implementation details. Additionally we propose a partially parallel version of the quilting technique.

A. Introduction

Texture synthesis is a classical image processing problem that finds its applications in virtual reality rendering (video games, animation movies, . . . ). Given an input texture image, it consists in producing output texture images that are both visually similar to and pixel-wise different from the input, and having possibly a larger size. One can separate texture synthesis algorithms into two categories, namely statistical constraint approaches and non-parametric patch-based methods, although "hybrid" algorithms have been proposed recently [] and the algorithm we propose in Chapter .

Statistical constraint texture synthesis algorithms model the texture based on statistical and/or perceptual considerations. They generally involve two different steps, one for analysis and one for synthesis. The analysis step consists in estimating a set of relevant statistics from the input texture image. The synthesis step computes an image that satisfies the statistical constraints estimated during the analysis step. Following the seminal paper of Heeger and Bergen [, ], several methods are based on statistics of wavelet coefficients or more involved multiscale image representations [, , ]. Another approach initially proposed by van Wijk [] consisting in randomizing the Fourier phase of an image has been extended to an exemplar-based synthesis method in [, ] suitable for micro-textures (textures that do not present salient geometric patterns and that are not constituted of individual discernible objects).

Non-parametric patch-based algorithms attempt at producing a new texture by arranging local neighborhoods of the input texture in a consistent way. The first methods of this kind are sequential algorithms that create a new texture one pixel at a time [, ]. These algorithms represented  a breakthrough in the field since they are able to reproduce macro-textures with specific geometric structures, see e.g. the IPOL paper []. Along with the method of Liang et al [] developed at the same period, image quilting [], the algorithm investigated in this paper, computes the new texture by arranging seamlessly full patches of the input texture. The main innovation of this method is the procedure to stitch a new patch in the sequentially built output texture to avoid discontinuities as much as possible. This is achieved by computing an optimal boundary cut between the patch and the synthesis area thanks to a linear programming optimization (see Section A. for details). Let us note that another solution of this "stitching step" is proposed in [?] using graphcuts. Many contributions have since improved the results of image quilting, at least regarding the computational cost, and we refer to the state of the art [] for a more complete survey of this category of texture synthesis algorithms. Let us also mention the recent paper [], based on [?, ], that discusses and atempts to solve some limitations of these approaches.

The plan of the paper is as follows. Section A. describe in detail all the steps of the algorithm and in particular gives a detailed mathematical description of the linear programing problem for the computation of the minimum error boundary cut. Section A. gives implementation details and discuss a parallelization of the Efros-Freeman algorithm which, to the best of our knowledge, is a contribution of this paper. Section A. presents numerous experiments of our implementation of the Efros-Freeman algorithm. This experimental section shows that this algorithm produces most of the time visually good results. However, a set of failure cases shows that the shortcomings of the Efros-Leung algorithm, that is, garbage growing and verbatim copy [, ], are also present in the Efros-Freeman results, but only at a larger scale.

A. Algorithm Description

In [] the authors propose a sequential patch-based algorithm to synthesize textures. I 0 and I s denote the input sample and the output texture respectively. The output image I s is constructed patch by patch in a raster scan order. The goal of each iteration is to fill a patch P old of I s that is only partially defined on a region called overlap region (see Figure A.). To do so a patch P in of I 0 that matches P old on the overlap region is randomly selected. An optimal boundary cut between P old and P in is then computed within the overlap region. This optimal boundary cut is used to construct the new patch P new by blending P old and P in along the cut.

The whole image quilting algorithm is described in Algorithm  and the remaining of this section will detail each step of the algorithm, namely the initialization, the patch search procedure to select P in , the computation of the minimum error boundary and the blending procedure along the boundary.

A.. Initialization

The first step of the algorithm is to initialize I s . For that a random patch P in of size w p × w p is taken from I 0 and placed at the top-left corner of I s .

A.. Patch Selection

Once the image has been initialized the algorithm synthesizes the remaining patches of I s sequentially in raster scan order. At each iteration one has to fill a patch P old of I s that is only defined on an overlap region of width w o . Note that there are three possible overlap regions: vertical overlap Compute minimum error boundary cut between P old and P in (see Algorithm )

:

Construct the patch P new by blending P old and P in along the boundary cut (see Equation (A.))

:

Replace P old with P new within I s : end for To select a patch P in of the input image I 0 one computes the square distance between the overlap region of the patch P old of I s and the corresponding regions of all the patches of I 0 . The minimal distance d min is determined and P in is randomly picked among all patches whose distance to P old is lower than (1 + ε)d min where ε is the tolerance parameter.

To conclude this section let us give a detailed expression of the distance used to compare patches. A patch of I 0 is represented by the position of its top-left corner (m, n) ∈ {0, . . . , M 0w p } × {0, . . . , N 0w p }. The squared distance image D contains at each position (m, n) the distance between P old and the patch from I 0 who's top-left corner is (m, n) according to some binary weight Q that equals one in the overlap region and zero otherwise. More precisely for all (m, n) ∈ {0, . . . , M 0w p } × {0, . . . , N 0w P }, one has

D(m, n) = wp-1 i,j=0 Q(i, j)(P (i, j) -I 0 (m + i, n + j)) 2 . (A.)
The computation of the squared distance image D is discussed in Section A.. and the whole patch selection procedure is summarized in Algorithm .



A.. Minimum Error Boundary Cut

This step of the algorithm is the main contribution of the Efros-Freeman algorithm [].

The Patch Search step (see Section A..) gives a patch P new of I 0 having coordinates (m, n) that is similar to the partially defined current patch P old in their overlap region. We recall that the overlap regions that arise in the raster scan order are of three types: vertical, horizontal, and L-shaped overlap. These three cases are illustrated in Figure A.. To get the final patch P one must combine the patches P old and P new . Denoting Q the binary weight for the overlap regions as in the previous section, then, for any binary image M such that 0 ≤ M ≤ Q, P can be defined as the combination

P new = M P old + (1 -M ) P in .
The main idea of Efros and Freeman [] is looking for a binary shape M where the transition between P old and P new along the boundary of the shape is minimal. For simplicity, and to be able to use linear programming, the authors do not allow for any shape, but only for the ones whose boundaries are simple forward paths from one end to the other of the overlap region. This results in two pieces of image that are sewn together along some general boundary path, hence the algorithm name "quilting". In the remaining of this section, we describe in details the computation of the minimum error boundary cut for the three overlap cases, starting with the vertical and horizontal cases.

Connected Paths and Boundary Error

In the following, we call a path of length ℓ ≥ 1 in {0, . . . , w p -1} 2 any ordered sequence γ

= (γ 0 , γ 1 , . . . , γ ℓ-1 ) of -connected pixels γ k ∈ {0, . . . , w p -1} 2 , that is, for all k ∈ {0, . . . , ℓ -2}, max(|γ 1 k+1 -γ 1 k |, |γ 2 k+1 -γ 2 k |) = 1
(for convenience, within this section the first and second coordinates of a pixel γ k are denoted by γ 1 k and γ 2 k ). Let us recall that we use matrix coordinates for the pixels. Hence a path γ of length ℓ is said to be vertical if for all k ∈ {0, . . . , ℓ -2}, γ 1 k+1γ 1 k = -1 (vertical paths are oriented from bottom to top), and horizontal if for all k ∈ {0, . . . , ℓ -2}, γ 2 k+1γ 2 k = -1 (horizontal paths are oriented from right to left).

Denoting by e(i, j) = (P new (i, j)-P old (i, j)) 2 the squared difference between the two patches P old and P new , the boundary error of a path γ of length ℓ is defined by

E(γ) = ℓ-1 k=0 e(γ k ).



Vertical Boundary Cuts

In this section we explain how the optimal boundary between P old and P new is defined and computed in the case of vertical overlap, that is when the overlap region is the rectangle {0, . . . , w p -1} × {0, . . . , w o -1}.

The optimal boundary is defined as the vertical path γ that minimizes the boundary error E(γ) while joining both ends of the overlap regions. More precisely, define the admissible vertical paths as

Γ v = {γ = (γ 0 , . . . , γ wp-1 ), ∀k, γ 1 k = w p -1 -k and γ 2 k ∈ {0, . . . , w o -1}} (one notices that paths of γ ∈ Γ v are indeed vertical since γ 1 k+1 -γ 1 k = -1).
Then an optimal boundary is defined as any solution of the optimization problem min γ∈Γv E(γ).

(A.)

Problem (A.) is solved using dynamic programming. This is possible because Problem (A.) verifies the principle of optimality: if (γ 0 , . . . , γ wp-1 ) is an optimal solution of Problem (A.), then, for all 0 ≤ r ≤ w p -1, the subpath (γ 0 , . . . , γ r-1 ) is an optimal vertical path to reach the r-th point γ r-1 starting from the bottom of the overlap region.

Let us now discuss in detail the dynamic programming method to solve Problem (A.). For all points (i, j) ∈ {0, . . . , w p -1} × {0, . . . , w o -1}, denote by Γ (i,j) v the set of vertical paths γ that start at the bottom side of the overlap region and end at the point (i, j), that is

Γ (i,j) v = {γ = (γ 0 , . . . , γ wp-1-i ), ∀k, γ 1 k = w p -1-k, γ 2 k ∈ {0, . . . , w o -1}, and γ wp-1-i = (i, j)}.
Now, for all points (i, j) ∈ {0, . . . , w p -1} × {0, . . . , w o -1}, define E v (i, j) as the minimal cumulative vertical error to reach (i, j) starting from the bottom side, that is,

E v (i, j) = min γ∈Γ (i,j) v E(γ).
Then, since

Γ v = j∈{0,...,wo-1} Γ (0,j) v , one has min γ∈Γv E(γ) = min j∈{0,...,wo-1} E v (0, j).
Now, remark that the last-but-one point of an optimal vertical path that ends at (i, j) is one of the (at most) three points below (i, j), namely the points (i + 1, j -1), (i + 1, j), (i + 1, j + 1) (there can be only two neighboring points if (i, j) is at the border, that is, if j = 0 or j = w o -1). Hence, for all i ∈ {0, . . . , w p -2}, one has

E v (i, j) = e(i, j) + min(E v (i + 1, j -1), E v (i + 1, j), E v (i + 1, j + 1)), (A.)
where by convention, if one of the index j -1 or j + 1 is not a valid index, the corresponding term E v (i + 1, j -1) or E v (i + 1, j + 1) is discarded from the minimum. Besides, for the last line i = w p -1, one simply has E v (w p -1, j) = e(w p -1, j) since the paths are only made of one pixel. Hence the dynamic programming procedure for solving Problem (A.) is to compute the costs E v (i, j) line by line from bottom to top using Equation (A.) recursively, and then search for the minimal value of the first line j * = arg min j E v (0, j). The full optimal path γ can then  be traced back starting from this coordinate (0, j * ) = (γ 1 wp-1 , γ 2 wp-1 ). More precisely, if (i, j) are the coordinates of the point γ wp-1-i of the optimal path γ, then its preceding point γ wp-1-(i+1) is given by

γ wp-1-(i+1) = arg min(E v (i + 1, j -1), E v (i + 1, j), E v (i + 1, j + 1)).
In practice one stores the matrix T v (i, j) = arg min(E v (i+1, j -1), E v (i+1, j), E v (i+1, j +1)) while computing the vertical cumulative error E v given by Equation (A.) in order to trace back the path without additional computation since for i = w p -2 to 0, γ i = T v (γ i+1 ).

Horizontal Boundary Cuts

The case of horizontal overlap is just the symmetric case of vertical overlap one. Still let us introduce the notation that will be necessary for dealing with the L-shaped overlap case. One defines the set of horizontal paths as

Γ h = {γ = (γ 0 , . . . , γ wp-1 ), ∀k, γ 2 k = w p -1 -k and γ 1 k ∈ {0, . . . , w o -1}},
and for all points (i, j) ∈ {0, . . . , w o -1} × {0, . . . , w p -1}, one defines

Γ (i,j) h = {γ = (γ 0 , . . . , γ wp-1-j ), ∀k, γ 2 k = w p -1-k, γ 1 k ∈ {0, . . . , w o -1}
, and γ wp-1-j = (i, j)} the set of paths that start at the right side of the overlap region and end at (i, j), as well as

E h (i, j) = min γ∈Γ (i,j) h E(γ)
the minimal cumulative horizontal error to reach (i, j) starting from the right side of the overlap region. Then E h (i, w p -1) = e(i, w p -1), and for all j ∈ {0, . . . , w p -2}, one has the recursive relation E h (i, j) = min(E h (i -1, j + 1), E h (i, j + 1), E h (i + 1, j + 1))

which enables to compute E h (i, j) for all (i, j) ∈ {0, . . . , w o -1} × {0, . . . , w p -1}. Since

min γ∈Γ h E(γ) = min i∈{0,...,wo-1} E h (i, 0)
the end point of the optimal horizontal path is (i * , 0) where i * = arg min i E h (i, 0) and the path can be traced back thanks to the matrix T h (i, j) = arg min(E h (i -1, j + 1), E h (i, j + 1), E h (i + 1, j + 1)).

L-shaped Boundary Cuts

The case of L-shaped overlap regions is slightly more complex than the vertical and horizontal cases since the geometry of the L-shape does not allow for a clear ordering of the pixels of the path. The original paper [] only mentions that, in the L-shaped overlap case, "the minimal paths meet in the middle and the overall minimum is chosen for the cut". We propose below a rigorous interpretation of this sentence.

A separating path will start at the bottom side of the L-shape (that is, the points (i, j) with i = w p -1 and j ∈ {0, . . . , w o -1}) and end at the right side of the L-shaped (that is, the points (i, j) with i ∈ {0, . . . , w o -1} and j = w p -1). To restrict the geometry of admissible paths (and allowing for dynamic programming), it is further required that an L-shaped path has to meet  at some diagonal point (i, i), i ∈ {0, . . . , w o -1}, that the first part from the bottom side to this diagonal point (i, i) is a vertical path, and that the remaining part from the diagonal point (i, i) to the right side is a reversed horizontal path. More formally, for all indexes i ∈ {0, . . . , w o -1}, one defines

Γ i L = {γ = (γ 0 , . . . , γ 2(wp-i)-1 ), (γ 0 , . . . , γ wp-i-1 ) ∈ Γ (i,i) v and (γ 2(wp-i)-1 , . . . , γ wp-i-1 ) ∈ Γ (i,i) h },
and for all γ ∈ Γ i L we denote by γ v and γ h its associated vertical and horizontal paths of

Γ (i,i) v and Γ (i,i) h
respectively. The set of the admissible L-shaped boundaries Γ L is then defined as the disjoint union of the sets Γ i L , that is,

Γ L = w 0 -1 i=0 Γ i L .
As before, we search for an optimal L-shaped boundary cut γ ∈ Γ L having minimal boundary error

E(γ) = ℓ(γ)-1 k=0 e(γ k ),
where ℓ(γ) is the length of the path γ (which is equal to 2

(w p -i) -1 if γ ∈ Γ i L )
. The optimal L-shaped boundary can be found in splitting the above sum into a vertical part and an horizontal part, since for all γ ∈ Γ i L one has

E(γ) = E(γ v ) + E(γ h ) -e(i, i).
Hence one has

min γ∈Γ L E(γ) = min i∈{0,...,wo-1} min γ∈Γ i L E(γ) = min i∈{0,...,wo-1} min γ∈Γ i L E(γ h ) + E(γ v ) -e(i, i) = min i∈{0,...,wo-1} min γ v ∈Γ (i,i) v E(γ v ) + min γ h ∈Γ (i,i) h E(γ h ) -e(i, i) = min i∈{0,...,wo-1} E v (i, i) + E h (i, i) -e(i, i),
where E v and E h are the vertical and horizontal cumulative errors defined in the previous sections. Hence,

min γ∈Γ L E(γ) = min i∈{0,...,wo-1} E v (i, i) + E h (i, i) -e(i, i).
The above equation enables us to determine the optimal position (i * , i * ) on the diagonal of the optimal path γ once the matrices E v and E h have been computed recursively. We can then trace back the vertical part γ v and the horizontal part γ h of the optimal path γ using the matrices T v and T h . Let us remark that the "overall minimum" evoked in [] must not be interpreted as the minimum of

E v (i, i) + E h (i, i) but the minimum of E v (i, i) + E h (i, i) -e(i, i).
without normalization. Hence the operators that are computed are respectively M 0 N 0 F for the forward transform and F -1 for the backward transform. Hence we have

V * W = F -1 (F(V * W )) = 1 M 0 N 0 F -1 (M 0 N 0 F(V )M 0 N 0 F(W )).
This means that after multiplying the two FFTW forward transforms and performing the backward inverse transform, one has to normalize in dividing by the size of the images.

Cross-correlation

Let us denote by Γ(V, W ) the cross-correlation between two images,

Γ(V, W )(m, n) = M 0 -1 k=0 N 0 -1 l=0 V (k, l)W (m + k, n + l).
Note also Ṽ the symmetric of V with respect to the origin, that is,

Ṽ (m, n) = V (-m, -n). Note that one has Γ(W, V )(m, n) = Γ(V, W )(-m, -n), that is, Γ(W, V ) = Γ(V, W
) and that the cross-correlation is simply a convolution between the symmetric version of the first image and the second image, that is,

Γ(V, W )(m, n) = Ṽ * W (m, n) = V * W (-m, -n).
Hence the computational cost for a cross-correlation image is the same as the one for a convolution product, that is three FFT calls. We need to compute the squared distance between a patch P old and all the patches of I 0 according to some binary weight Q that represents the overlap region (that is a horizontal, vertical, or L-shaped mask). More precisely for all (m, n) ∈ {0, . . . , M 0w p } × {0, . . . , N 0w p }, we want to compute D(m, n) = wp-1 i,j=0 Q(i, j)(P old (i, j) -I 0 (m + i, n + j)) 2 .

Let us denote by P ext and Q ext the extensions of P old and Q into images of size M 0 × N 0 by filling the domain with 0-valued pixels, that is Q(i, j)P old (i, j) 2 -2Γ(Q ext P ext , I 0 )(m, n) + Γ(Q ext , I 2 0 )(m, n),



where the multiplications Q ext P ext and I 2 0 are componentwise. Within our procedure, we will compare a patch P old that is taken from the output image that is only partially defined, and the binary weight Q representing the overlap region. Hence, if all the undefined pixels of the output patch P old are set to 0 (by initializing I s to 0), one has ∀(i, j) ∈ {0, w p -1} 2 , Q(i, j)P old (i, j) = P old (i, j).

Consequently,

wp-1 i,j=0 Q(i, j)P old (i, j) 2 = wp-1 i,j=0 P old (i, j) 2 = P old 2 2 and Q ext P ext = P ext , and thus Γ(Q ext P ext , I 0 ) = Γ(P ext , I 0 ). Hence, the binary mask Q is only influent when computing Γ(Q ext , I 2 0 ), a computation that is done just once at the beginning of the procedure since it does not depend on P old .

In the end, the image D is simply given by

D = P old 2 2 -2Γ(P ext , I 0 ) + Γ(Q ext , I 2 0 ) (A.)
where the last image Γ(Q ext , I 2 0 ) is computed once before running the algorithm and stored in memory. Hence the cost for computing D is only three FFT calls for the computation of the crosscorrelation Γ(P ext , I 0 ), and is even reduced to two FFT calls by storing the DFT of I 0 .

Note that the image D is defined for all (m, n) ∈ {0, . . . , M 0 -1} × {0, . . . , N 0 -1}, but the squared distances for (m, n)-coordinates outside {0, . . . , M 0 -w p }×{0, . . . , N 0 -w p } correspond to patches of size w p × w p that are defined by periodic boundary conditions and thus those patches must be discarded when searching for the minimal distance.

Minimum distance

Theoretically the results of the distance computation using the sum square differences or the FFT are equal. But in practice, the computation using the FFT is subject to rounding errors. Especially if there is a patch such that D(m, n) = 0, when adding positive and negative terms in Equation (A.) the result might be negative due to numerical errors. This rounding error problem leads to a negative minimal squared distance d min , which leads to errors in the patch search algorithm (Algorithm ) since then none of patches satisfy the condition D(m, n) < (1 + ε)d min . To avoid this, all squared distances D(m, n) smaller than 1 are set to D(m, n) = 1 (the value 1 is smallest distance-value greater than 0).

A.. Parallelization

Although the Efros-Freeman algorithm was presented as a purely sequential algorithm, we found out that the procedure can be partially parallelized, and thus the algorithm significantly accelerated when running on multi-core platforms.

In general, non-parametric patch based methods can not be completely parallelized because the sequential assignment of patches is strongly dependent of the previous synthesis steps. Nevertheless, in the case of the quilting method, a new patch assignment is only correlated to some previous synthesis steps and not all of them. To illustrate this let us consider the first few steps of the algorithm. The first step of the algorithm adds the seed patch and the second step adds a second patch on the right. When adding the third patch, the first patch of the second row can be created at the same time since it only depends on the one that is on top of it and on its top-right. More generally the synthesis of a new row of patches can be started as soon as two patches of the previous row have been synthesized. This procedure is illustrated in Let us discuss the acceleration induced by the parallelization. Suppose we synthesize an image I s made of N r × N c patches where N r is the number patch rows and N c is the number of patch columns. Without the parallelization the number of iterations N it is N it = N r N c . When parallelizing the number of iterations N it is reduced to 2(N r -1) + N c . This parallelization is especially effective if one has as many processors available as patches to synthesize simultaneously in one iteration. The maximal number of patches to synthesize simultaneously is equal to min (N r , Nc+1

2

), where ⌊x⌋ = max(n ∈ N, n ≤ x), since within an iteration there is at most one new patch per row and one new patch per pair of column.

The parallel version of the quilting algorithm is summarized in Algorithm . Let us notice that for this parallelization to be valid the overlap size must satisfy w o ≤ w p /2 which is always the case in practice. end for : end for 

Algorithm : Parallelization

A. Experiments

In this section texture synthesis results are shown using the algorithm described previously []. Two apects of this method are emphasized. On the one hand the visual quality result and on the other hand the tendency of the method to verbatim copy parts of the input sample. For this we represent each result as shown in Figure A.: the input sample I 0 , the color map, the synthesized texture I s and the synthesis map. The color map is an image of size equal to the size of I 0 and each of its pixels is assigned a different color. This allows to visualize each position of I 0 with the corresponding color of the color map. The synthesis map allows to identify from where comes every patch copied in I 0 and quilted in I s . This allows to easily identify verbatim copies. Each pixel p ′ i in I s is assigned a pixel p j from I 0 . The synthesis map at position p ′ i is mapped to the value of the color map at position p j .

Let us discuss the influence of the parameters. There are three of them: the patch size w p , the overlap size w o and the tolerance parameter ε. The overlap size is expressed as the proportion respect to the patch size. That is w o = 0.25 implicitly means w o = 0.25w p .

In Figures A. The last failure case is the "growing garbage" drawback. This is when the method is stuck in a region of the input sample repeating it on a large part of the output. This is due to the local aspect of the algorithm and it is more noticeable when the texture examples are not stationary. This is also shown in Figure A. in the last row. It is important to notice that the two last failure cases are more noticeable when the output size -input size ratio is higher than two.

Influence of the patch size



smaller patches this effect is reduced, since more "similar" patches are available in patch search step and then the set of candidate patches is larger allowing more variation.

Influence of the overlap size

The overlap size is an important parameter in this method. To analyze its influence w p = 40 and ε = 0.1 are fixed and the results for w o ∈ {0.10, 0.25, 0.5} are compared. The general conclusion is that for larger overlap sizes the transition from one patch to the other is satisfying on the cost of growing garbage for some texture sample as can be seen in the two first examples in

Influence of the tolerance parameter

This last parameter is directly related to the verbatim copy effect. For this analysis w p = 40 and w o = 0.25 are fixed and the results for ε ∈ {0.05, 0.1, 0.3, 0.5, 0.7} are compared. Increasing ε implies having more candidate patches for the synthesis. This allows more variation when choosing a patch in the image and this is directly seen in the synthesis maps of the examples in Figure A. where for ε ∈ {0.5, 0.7} the patches are taken more "randomly" from the input sample thus reducing the size of the verbatim copy regions. For some texture examples the visual quality of the result can decrease. On the other hand, as expected, low values of ε leads to very large verbatim copy areas.

A. Conclusion

In this paper we analyzed in detail Efros and Freeman's texture synthesis algorithm []. Extensive numerical experiments have been proposed to illustrate the performance of the method as well as the influence of its parameters. We conclude from these experiments that in general, for the correct set of parameters, the visual results are satisfying to the cost of verbatim copying large parts of the input textures that might be visually disturbing. This is a common drawback of Efros and Leung's method []. Another issue that arises in [] is the garbage growing effect. This is especially apparent when the input texture is not stationary. We also noticed that it is related to the raster scan order used to synthesize the image, which propagates the errors. All these drawbacks are in general more apparent when synthesizing an image significantly larger than the input. Hence, due to the raster scan order, the quilting algorithm is not suitable to generate very large texture images since the quality tends to decrease with the distance to the top left corner of the image. We provide with this analysis a strategy to partially parallelize the method who was initially introduced as essentially sequential. This allows a significant speed up when running with multicore processors. Résumé : Cette thèse s'attaque au problème de la synthèse de texture par l'exemple en utilisant des modèles stochastiques locaux de patchs pour générer de nouvelles images.



La synthèse de texture par l'exemple a pour but de générer à partir d'un échantillon de texture de nouvelles images qui sont perceptuellement équivalentes à celle de départ. Les méthodes peuvent se regrouper en deux catégories: les méthodes paramétriques et les non paramétriques à base de patchs. Le premier groupe a pour but de caractériser une image de texture à partir d'un ensemble de statistiques qui définissent un processus stochastique sous-jacent. Les résultats visuels de ces méthodes sont satisfaisants, mais seulement pour un groupe réduit de types de texture. La synthèse pour des images de textures ayant des structures très contrastées peut échouer. La deuxième catégorie d'algorithme découpe, puis recolle de manière consistente des voisinages locaux de l'image de départ pour générer de nouvelles configurations plausibles de ces voisinages (ou patchs). Les résultats visuels de ces méthodes sont impressionnants. Néanmoins, on observe souvent des répétitions verbatim de grandes parties de l'image d'entrée qui du coup peuvent être reproduites plusieurs fois. De plus, ces algorithmes peuvent diverger, reproduisant de fac ¸on itérative une partie de l'image de l'entrée en négligeant le reste.

La première partie de cette thèse présente une approche combinant des idées des deux catégories de méthodes, sous le nom de synthèse localement Gaussienne. On préserve dans cette nouvelle méthode les aspects positifs de chaque approche: la capacité d'innover des méthodes paramétriques, et la capacité de générer des textures fortement struc-turées des méthodes non paramétriques à base de patchs. Pour ce faire, on construit un modèle Gaussien multidimensionnel des auto-similarités d'une image de texture. Ainsi, on obtient des résultats qui sont visuellement supérieurs à ceux obtenus avec les méthodes paramétriques et qui sont comparables à ceux obtenus avec les méthodes nonparamétriques à base de patchs tout en utilisant une paramétrization locale de l'image. La thèse s'attache aussi à résoudre une autre difficulté des méthodes à base de patchs: le choix de la taille du patch. Afin de réduire significativement cette dépendance, on propose une extension multi échelle de la méthode.

Les méthodes à bases de patchs supposent une étape de recollement. En effet, les patchs de l'image synthétisée se superposent entre eux, il faut donc gérer le recollement dans ces zones. La première approche qu'on a considérée consiste à prendre en compte cette contrainte de superposition dans la modélisation des patchs. Les expériences montrent que cela est satisfaisant pour des images de textures périodiques ou pseudo-périodiques et qu'en conséquence l'étape de recollement peut être supprimée pour ces textures. Cependant, pour des images de textures plus complexes ce n'est pas le cas, ce qui nous a menée à suggérer une nouvelle méthode de recollement inspirée du transport optimal.

Cette thèse conclut avec une étude complète de l'état de l'art en génération d'images de textures naturelles. L'étude que nous présentons montre que, malgré les progrès considérables des méthodes de synthèse à base d'exemples proposées dans la vaste littérature, et même en les combinant astucieusement, celles-ci sont encore incapables d'émuler des textures complexes et non stationnaires. Exemplar based texture synthesis is the process of generating, from an input texture sample, new texture images that are perceptually equivalent to the input. There are roughly two main categories of algorithms: the statistics-based methods and the nonparametric patch-based methods. The statisticsbased methods aim at characterizing a given texture sample by estimating a set of statistics which will define an underlying stochastic process. The new images will then be samples of this stochastic process, i.e. they will have the same statistics as the input sample. The question here is what would be the appropriate set of statistics to yield a correct synthesis for the wide variety of texture images? The results of statistical methods are satisfying but only on a small group of textures, and often fail when important structures are visible in the input. The non-parametric patch-based methods reorganize local neighborhoods from the input sample in a consistent way to create new texture images. These methods return impressive visual results. Nevertheless, they often yield verbatim copies of large parts of the input sample. Furthermore, they can diverge, starting to reproduce iteratively one part of the input sample and neglecting the rest of it, thus growing what experts call "garbage". In this thesis we propose a technique combining ideas from the statistical based methods and from the non-parametric patch-based methods. We call it the locally Gaussian method. The method keeps the positive aspects of both categories: the innovation capacity of the parametric methods and the ability to synthesize highly structured textures of the non-parametric methods. To this aim, the selfsimilarities of a given input texture are modeled with conditional multivariate Gaussian distributions in the patch space. A new image is generated patchwise, where for each given patch a multivariate Gaussian model is inferred from its nearest neighbors in the patch space of the input sample, and hereafter sampled from this model. The synthesized textures are therefore everywhere different from the original. In general, the results obtained are visually superior to those obtained with statistical based methods, which is explainable as we use a local parametric model instead of a global one. On the other hand, our results are comparable to the visual results obtained with the non-parametric patch-based methods. This dissertation addresses another weakness of patch-based methods. They are strongly dependent on the patch size, which has to be decided manually. It is therefore crucial to fix a correct patch size for each synthesis. Since texture images have, in general, details at different scales, we extend the method to a multiscale approach which reduces the strong dependency of the method on the patch size. Patch based methods involve a stitching step. Indeed, the patches used for the synthesis process overlap each other. This overlap must be taken into account to avoid any transition artifact from patch to patch. Our first attempt to deal with it was to consider directly the overlap constraints in the local parametric model. The experiments show that for periodic and pseudo-periodic textures, considering these constraints in the parametrization is enough to avoid the stitching step. Nevertheless, for more complex textures it is not sufficient. This leads us to suggest a new stitching technique inspired by optimal transport and midway histogram equalization. This thesis ends with an extensive analysis of the generation of several natural textures. This study shows that, in spite of remarkable progress for local textures, the methods proposed in the extensive literature of exemplar based texture synthesis are still limited when dealing with complex and nonstationary textures.

FIGURE 1 -

 1 FIGURE 1 -On montre trois it érations diff érentes de la synth èse d'une texture. A chaque it ération un patch est synth étis é. Il est repr ésent é en rose à chaque it ération. De gauche à droite on observe les trois cas d'intersection : vertical, horizontal et en forme de L.

Chapitre  :

  FIGURE 2 -Influence des param étres. L'image du coin sup érieur gauche et la texture d'entr ée. Les images restantes sont des synth éses en variant la taille du patch n et m la taille de l'ensemble U . De gauche à droite n = 10, 20, 30. Du haut vers le bas m = 10, 20, 30. Tous ces r ésultats correspondent à une intersection o de taille n/2.

FIGURE 3 -

 3 FIGURE 3 -Comparaison des modèles de patchs. De gauche à droite : texture d'entr ée, synth èse LG, synth èse CLG et synth èse RCLG. Aucune étape de recollement a ét é appliqu ée pour les trois mod èles montr ès. Les param ètres pour ces synth èses sont n = 40, m = 30 et o = 0.5n.

FIGURE 4 -

 4 FIGURE 4 -Comparaison de la synthèse localement Gaussienne (LG) à la synthèse multi-échelle (M SLG). De gauche à droite : texture d'entr ée, synth èse LG et synth èse M SLG. Les param ètres utilis és pour ces synth èses sont n = 20, m = 50 et o = 0.5n. Pour le multi-échelle le nombre d' échelles K utilis ées est K = 3. Pour les deux mod èles on a utilis é le m ême patch d'initialisation.

FIGURE 6 -

 6 FIGURE 6 -Illustration des signaux de dimension un à recoller. Les patchs à recoller sont repr ésent és par I 1 et I 2 et les bords de recollement sont repr ésent és par s 1 et s 2 .

  FIGURE 8 -Texture de bois. L'image originale est de taille 7087 × 23622 pixels. Pour les exp ériences on a travaill é sur des images r éduites d'un facteur 16.

  FIGURE 10 -R ésultats de synth èse pour diff érentes "grande" textures.

Figure 1

 1 Figure 1.1 -Three different iterations of the synthesis process are shown. At each iteration a patch is being synthesized. This patch is represented by the pink square in the three iterations shown. From left to right the three overlap cases are represented: vertical, horizontal and L-shape.

  Figure 1.2 -Texture synthesis result for the left top corner texture image. We show the results obtained for different values of m (the number of similar patches) and n (the side patch size). From left to right n = 10, 20, 0. From top to bottom, the number of nearest neighbours is m = 10, 20, 30. All the results are obtained for an overlap of a half patch size o = n/2.

Figure 1 . 3 -

 13 Figure 1.3 -Patch model comparison. From left to right: texture sample, synthesis result using LG, CLG and RCLG. No quilting technique was applied to stitch together the simulated patches for the three presented models. The parameters used for all examples are n = 40, m = 30 and o = 0.5n.

Figure 1 . 4 -

 14 Figure 1.4 -Comparison of LG (Locally Gaussian) to MSLG (Multiscale Locally Gussian ). From left to right: texture sample, synthesis result using LG, and MSLG. The parameters used for these examples are n = 20, m = 50 and o = 0.5n. For the multiscale synthesis results the number of scales K used is K = 3.For both methods we used the same seed patch for each texture sample.

Figure 1

 1 Figure 1.5 -Teaser. From left to right: input, input blended with Poisson editing [47], input blended with quilting [15], input blended with midway blending.

Figure 1

 1 Figure 1.6 -Illustration of the blending edges. The two images to blend are represented by I 1 and I 2 and the edges on which the midway blending is applied are s 1 and s 2 .

Figure 1

 1 Figure 1.7 -Some challenging examples blending patches which do not match. Comparison with two other stitching methods. From left to right: input, Poisson editing [47] and the midway blending.

  Figure . and its crops in Figure .. We follow by showing results of a combination of the considered texture synthesis methods, namely:

Figure 1 Figure 1

 11 Figure 1.8 -Wood texture 1. The original input provided was of size 7087 × 23622 pixels. For the experiments the image was reduced by a zoom out of factor 16.

  Figure . we show three examples.

Figure 1

 1 Figure 1.10 -Synthesis results of different natural large textures.

u

  Ω → R: input texture image defined on the discrete domain Ω = I M × I N of size M × N where I c denotes the discrete interval [0, . . . , c -1] w Ω r → R: output texture image defined on the discrete domain Ω r = I rM × I rN of size rM × rN d is equal to 1 if u is grayscale and to 3 for color images r ratio (size of output image w)/(size of input image u) n side patch size m number of nearest neighbours used to learn the parameters of the Gaussian distribution o overlap size p (x,y) u square patch of size n × n from an image u of size M × N at position (x, y), p (x,y) u

Algorithm  :

  Naive locally Gaussian texture synthesis Input: input texture sample u, side patch size n, number of nearest neighbours m, overlap size o Output: reconstructed texture w : for x = 1 : n(1o) : M do : for y = 1 : n(1o) : N do :

  The size n of the patches and the size of the neighborhood m are user-specified values. The results of this first model are shown in Figure .. Observe that replacing the patches of the input texture by simulated ones (with the corresponding Gaussian models) achieves a faithful random variant of the original image, for reasonable values of n and m. For large values of m the model

  input m = 10, n = 20 m = 10, n = 30 m = 20, n = 10 m = 20, n = 20 m = 20, n = 30 m = 30, n = 10 m = 30, n = 20 m = 30, n = 30

Figure 2 . 1 -

 21 Figure 2.1 -Visual evaluation of the validity of a Gaussian patch model. This experiment performs a Gaussian substitution for each patch of the left top corner texture image (Algorithm 1). From left to right the patch size is n = 10, 20, 30. From top to bottom the number of nearest neighbours is m = 10, 20, 30.

  Figure ., two texture examples are presented.

Figure 2

 2 Figure 2.2 -Behaviour of the mean standard deviation of the patches as functions of the patch size n and the neighbourhood size m. For each texture example two patches have been selected (left column). For each of them the mean standard deviation is computed for several values of n ∈ {10, 20, 40, 80} and m ∈ {5, 10, 30, 50, 500, 1000}. In the input image patch 1 is represented by a solid line and patch 2 by a dashed line.

  sample of this model. The overlap size o is fixed and for all the presented results o is half the size n of the patch.

Figure 2

 2 Figure 2.3 -Three different iterations of the synthesis process are shown. At each iteration a patch is being synthesized. This patch is represented by the pink square in the three iterations shown. From left to right the three overlap cases are represented: vertical, horizontal and L-shape.

  input m = 10, n = 20 m = 10, n = 30 m = 20, n = 10 m = 20, n = 20 m = 20, n = 30 m = 30, n = 10 m = 30, n = 20 m = 30, n = 30

Figure 2 . 4 -

 24 Figure 2.4 -Texture synthesis result for the left top corner texture image. We show the results obtained for different values of m (the number of similar patches) and n (the side patch size). From left to right n = 10, 20, 0. From top to bottom, the number of nearest neighbours is m = 10, 20, 30. All the results are obtained for an overlap of a half patch size o = n/2.

Figure 2 . 5 -

 25 Figure 2.5 -Results representation. From left to right: texture sample, position map, synthesized image and synthesis map. The synthesis map shows for each synthesized patch its initial position in the texture sample. It allows then to identify exactly the verbatim copy regions (they correspond to continuous color areas of the map).

Figure 2 . 6 -

 26 Figure 2.6 -Comparison with various other texture synthesis methods. From left to right: texture sample, synthesis results of [20], [6], [50], [15] and our algorithm the locally Gaussian synthesis method (LG). For all the examples m = 10, n = 30 except for the third row where n = 10.



  Conclusions are presented in Section .. u Ω → R: input texture image defined on the discrete domain Ω = I M × I N of size M × N where I c denotes the discrete interval [0, . . . , c -1] w Ω r → R: output texture image defined on the discrete domain Ω r = I rM × I rN of size rM × rN d is equal to 1 if u is grayscale and to 3 for color images r ratio (size of output image w)/(size of input image u) n side patch size m number of nearest neighbours used to learn the parameters of the Gaussian distribution o overlap size K number of scales (maximum factor of zoom out is K

  (x,y) w is a vector of size k × 1 with the values of p (x,y) w on the overlap area and Sp (x,y) w is a vector of size  (dn 2k) × 1 with the other components of p (x,y) w

Figure 3

 3 Figure 3.1 -Three different iterations of the synthesis process are shown. At each iteration a patch is being synthesized. This patch is represented by the pink square in the three iterations shown. From left to right the three overlap cases are represented: vertical, horizontal and L-shape.

Figure 3

 3 Figure 3.2 -Comparison of LG to the quilting method [15]. These results are done taking equal parameters and seed patch to initialize the synthesized imaged. The parameters used are n = 40, o = n/2 and m = 30 (for the LG model).

  (x,y) w live in all R dn 2 , and  this ensures the existence of the conditional multivariate Gaussian distribution sought. The new covariance matrix is denoted by Γ = Σ + σ 2 I dn 2 .If we are in the case where the Gaussian vectors subspace p (x,y) w

  to conserve the low resolution synthesis from the previous scale. The high frequency patch p (x,y) Hw k is sampled form N (µ H , Σ H ) and then added to p (x,y)

  Figure ., general results of the multiscale method are shown with success and failure cases. In Figure . several texture synthesis methods are compared. In Figure . the innovation capacity of our method is compared to [] using coordinate maps. In Figure . the patch models introduced in Section . and . are compared. Finally, the influence of the parameters is discussed in Figure .. There are four of them: the patch size n, the number of neighbours m, the overlap size o and the number of scales K used in the multiscale approach. This is illustrated with two texture examples.

  Figure 3.3 -Synthesis results of the multiscale algorithm. The small images correspond to the texture samples and the big ones are the corresponding synthesis results. The parameters used for all examples are n = 20, m = 20 and K = 2.

Figure 3 . 4 -

 34 Figure 3.4 -Patch model comparison. From left to right: texture sample, synthesis result using LG, CLG and RCLG. No quilting technique was applied to stitch together the simulated patches for the three presented models. The parameters used for all examples are n = 40, m = 30 and o = 0.5n.

Figure 3 . 5 -

 35 Figure 3.5 -Comparison of LG to MSLG. From left to right: texture sample, synthesis result using LG, and MSLG. The parameters used for these examples are n = 20, m = 50 and o = 0.5n. For the multiscale synthesis results the number of scales K used is K = 3. For both methods we used the same seed patch for each texture sample.

  Figure . and . the influence of the parameters is illustrated. There are four of them: the patch size n ∈ {10, 20, 40}, the number of neighbours m ∈ {10, 30, 50}, the overlap size o ∈ {0.25, 0.5} × n and the number of scales used K ∈ {1, 2, 3}.

  This parameter corresponds to the number of patches used to estimate the patches' Gaussian distribution. As has been discussed in Section . the value of m in general should not be too small (> 5) or too large (< 50) to avoid patches of variance null or too big. These values are not general for all textures. The choice of m is linked to the amount of self-similarity in the image. Thus when m is too large the Gaussian sampling will blur up the image. In the experience of Figure . one can see that when m = 150 this leads to a texture which is too regular compared to the input sample. Otherwise since the sample texture in Figure . has many self similar patches the value of m can be large enough. The choice of m is a compromise between a copy-paste strategy and the risk of a blurry texture reconstruction.

Figure 3 . 6 -Figure 3

 363 Figure 3.6 -Comparison to several texture synthesis algorithms. From left to right: input sample, MSLG, MSLG combined to Portilla and Simoncelli [50], Portilla and Simoncelli[50], Tartavel et al. [59] (source: http://perso.telecom-paristech.fr/ ~tartavel/research/jmiv14.php) and Efros and Freeman [15]. For MSLG the parameters used are n = 20, m = 20 and K = 2. For [50] four scales and orientations were used. For [59] the patch size is 12 and the number of scales 3. For [15] the patch size used is 20.

Figure 3 . 9 -

 39 Figure 3.9 -Positive effect of a multiscale procedure. For each row from left to right: input sample, three synthesis results for K = 1, 2, 3 scales. For the three texture examples from top to bottom the parameters used were (n = 10, m = 20), (n = 20, m = 20), (n = 20, m = 30). In the three examples the global arrangement for K = 1 is not respected while for K = 2 and K = 3 are correctly synthesized.

Figure 4

 4 Figure 4.1 -Teaser. From left to right: input, input blended with Poisson editing [47], input blended with quilting [15], input blended with midway blending.

Figure 4

 4 Figure 4.2 -Illustration of the blending edges. The two images to blend are represented by I 1 and I 2 and the edges on which the midway blending is applied are s 1 and s 2 .

Figure 4

 4 Figure 4.3 -Illustration of the propagation band of width W . Left: image I 1 . The red rectangle marks out the propagation band of I 1 . The propagation direction is from bottom to top from A to B 1 . The edge in A represents the blending edge of I 1 . Center: image I 2 . The green rectangle marks out the propagation band of I 2 . The propagation direction is from top to bottom from A to B 2 . The edge in A represents the blending edge of I 2 . Right: output image Ĩ where the upper part is the corresponding portion in I 1 and the lower part is the corresponding portion of I 2 . The blue edge is the registered edge that matches perfectly after applying the midway deformations.

  Figure 4.4 -Midway blending of two synthetic images using the linear propagation method. From left to right: image I 1 , image I 2 , deformed image Ĩ1 , deformed image Ĩ2 , image I composed of the upper part of I 1 (above the red line) and the lower part of I 2 (under the red line) and Ĩ the blended version of I. Images Ĩ2 and Ĩ show the normalization artifact indicated within the red circles when propagating linearly the midway deformations.

  Figure .). For the propagation step each image is considered independently. Instead of working with the histograms we work with the inverse of the cumulative histograms. Thus for I 1 we consider H A and H B 1 the cumulative histograms of the edges in A and B 1 shown in Figure .. The edge in A is the blending edge after the midway deformation. The inverse cumulative histograms of the signals contained in the propagation band are a transition between H -1

Figure 4 . 5 -

 45 Figure 4.5 -Simple illustrative example 1. Left image: signals s 1 and s 2 to register. Middle image: cumulative histograms H 1 and H 2 and the midway histogram H. Right image: signals s 1 s 2 and the midway signal s. In this example one can observe that indeed s 1 and s 2 are translated to a signal that is halfway between them.

  Figure . where once again the behaviour of the midway deformation is not intuitive. In Figures . and . we consider two signals s 1 and s 2 having different mass. Once again one can observe that if the values of the flat zones are greater than zero then the deformed signals have and additional flat zone yielding unnatural merging results. In Figures . -. five blending results are shown for synthetic images. The propagation step is done using the progressive midway approach. The results obtained are coherent with the examples shown in the previous Figures . -..

Figure 4 . 7 -

 47 Figure 4.7 -Simple illustrative example 3. Left image: signals s 1 and s 2 to register. Middle image: cumulative histograms H 1 and H 2 and the midway histogram H. Right image: signals s 1 s 2 and the midway signal s.In this example one can observe that once again the deformation is not the expected translation. The signal is "split" in two and has an additional mode.

Figure 4 . 8 -Figure 4 Figure 4 Figure 4 Figure 4

 484444 Figure 4.8 -Simple illustrative example 5. Left image: signals s 1 and s 2 to register. Middle image: cumulative histograms H 1 and H 2 and the midway histogram H. Right image: signals s 1 s 2 and the midway signal s. This is the case of two signals with different mass. After renormalizing the midway signal is a translation of s 1 and s 2 to a half way signal.

Figure 4 .

 4 Figure 4.14 -This example illustrates a midway blending of two images I 1 and I 2 of different mass for a general case. The propagation of the deformation on each image is done using the progressive midway method on a band of 5 pixels of width. From left to right: inputs to blend, the deformed signals on I 1 and I 2 respectively, the blended output. Depending on the size of each flat regions and its values the resulting blended images can be as expected or not giving results similar to the ones obtained previously.

m 1 +m 2 2 :

 2 return s mid Algorithm : Midway blending  Input: Images to blend I 1 and I 2 , band width W , horizontal blending edge x 0  Output: Blended images Ĩ

Figure 4 .

 4 Figure4.15 -Illustration of the linear propagation of the deformation φ on a band of width W . This figure shows how the three points in A are modified after applying the deformation ψ(λ, t) yielding the three points in B. For the orange curve the point (x(t), y(t)) moves (x(φ(t)), y(φ(t))). For the blue curve the point (x λ (t), y λ (t)) moves (x λ (φ(λ, t)), y λ (φ(λ, t))). For the gray curve the point does not move.

Algorithm  :

  Generalized midway blending algorithm  Input: Images to blend I 1 and I 2 , band width W , reference position x 0 of the blending curve  Output: Blended images Ĩ : if edge-type = horizontal then :

FigureFigure 4

 4 Figure 4.16 -Different colour representations used to stitch two texture patches with the midway approach. From left to right: input, output using the RGB representation (notice the appearance of the colour artifacts along the stitching edge), output using the DCT representation and output using the PCA representation.

Algorithm  :

 : Texture synthesis  Input: Sample texture I, patch size n, overlap size o  Output: Synthesized texture Ĩ : Initialize Ĩ with a seed patch (randomly taken from I) : for each patch P A in Ĩ do :

  Figure .) or may present color artifacts (Figure .). For example in Figure . in the three fist examples one can clearly notice the quilting path. In Figure . the first example is an evident case that shows the advantages of the midway blending algorithm. This method tends to extend the structures when it can. Figure . shows failure examples.

Figure 4 .

 4 Figure 4.17 -Some challenging examples blending patches which do not match. Comparison with two other stitching methods. From left to right: input, Poisson blending [47], quilting method [15], midway blending.

Figure 4 .

 4 Figure 4.18 -Some challenging examples blending patches which do not match. Comparison with two other stitching methods. From left to right: input, Poisson blending [47], quilting method [15], midway blending.

Figure 4 .

 4 Figure 4.19 -Some challenging examples blending patches which do not match. Comparison with two other stitching methods. From left to right: input, Poisson blending [47], quilting method [15], midway blending.

Figure 4 .

 4 Figure 4.20 -Some nearly impossible examples blendig patches which do not match. Comparison with two other stitching methods. From left to right: input, Poisson blending [47], quilting method [15], midway blending.

Figure 4 .

 4 Figure 4.21 -Synthesis results 1. From left to right: input, texture synthesis results using midway blending, quilting algorithm [15]. For each example the patch size is n = 40, the same seed patch was taken to initialize the synthesized image, the overlap area is of size o = 0.45n and the band width is W = 10.

Figure 4 .

 4 Figure 4.22 -Synthesis results 2. From left to right: input, texture synthesis results using midway blending, quilting algorithm [15]. For each example the patch size is n = 40, the same seed patch was taken to initialize the synthesized image, the overlap area is of size o = 0.45n and the band width is W = 10.

Figure 4 .

 4 Figure 4.23 -Synthesis results 3. From left to right: input, texture synthesis results using midway blending, quilting algorithm [15]. For each example the patch size is n = 40, the same seed patch was taken to initialize the synthesized image, the overlap area is of size o = 0.45n and the band width is W = 10.



Figure 4 .

 4 Figure 4.24 -Synthesis results 4. From left to right: input, texture synthesis results using midway blending, quilting algorithm [15]. For each example the patch size is n = 40, the same seed patch was taken to initialize the synthesized image, the overlap area is of size o = 0.45n and the band width is W = 10.

Figure 4 .

 4 Figure 4.25 -Synthesis results 5. From left to right: input, texture synthesis results using midway blending, quilting algorithm [15]. For each example the patch size is n = 40, the same seed patch was taken to initialize the synthesized image, the overlap area is of size o = 0.45n and the band width is W = 10.

  , ., ., ., . different crops  of the inputs in Figures ., ., .

  , . and . which are crops of the inputs in Figures ., ., ., ., . and ..

  , ., ., ., ., ., ., .). This is clearly not the case in the smaller academic textures. Let us consider for instance the crops in Figure . which are very small portions (500 × 500 pixels) of the inputs in Figure . whose original size is of 7087 × 23633 pixels. The four crops

Figure 5

 5 Figure 5.1 -Wood texture 1. The original input provided was of size 7087 × 23622 pixels. For the experiments the image was reduced by a zoom out of factor 16.

Figure 5

 5 Figure 5.2 -Wood texture 2. The original input provided was of size 3500 × 13500 pixels. For the experiments the image was reduced by a zoom out of factor 16.

Figure 5

 5 Figure 5.3 -Wood texture 3. The original input provided was of size 4000 × 14000 pixels. For the experiments the image was reduced by a zoom out of factor 16.

Figure 5

 5 Figure 5.4 -Wood texture 4. The original input provided was of size 3400 × 14600 pixels. For the experiments the image was reduced by a zoom out of factor 16.

Figure 5

 5 Figure 5.5 -Stone textures 1, 2 and 3. The original inputs provided were of size 23622 × 23622, 6903 × 9983 and 21260 × 21260 pixels (left to right). For the experiments the images were reduced by a zoom out of factor 32, 8 and 32.

Figure 5

 5 Figure 5.6 -Metal texture 1 and 2. The original inputs provided were of size 23622 × 23622 and 6000 × 12390 pixels (left to right). For the experiments the images were reduced by a zoom out of factor 32 and 8 (left to right).

Figure 5

 5 Figure 5.7 -Four crops of different parts of Figure 5.1 at its original resolution. The cropped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

  Figure 5.8 -Four crops of different parts of Figure 5.2 at its original resolution. The cropped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

Figure 5

 5 Figure 5.9 -Four crops of different parts of Figure 5.3 at its original resolution. The cropped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

  Figure . four synthesis examples. The two first synthesis (top row in Figure .) show outstanding results. These microtextures are indeed well represented by their Fourier modulus. However this is not the case for the last two examples (bottom row in Figure .).

Figure 5 .

 5 Figure 5.10 -Four crops of different parts of Figure 5.4 at its original resolution. The cropped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

Figure 5 .

 5 Figure 5.11 -Four crops of different parts of Figure 5.5 at its original resolution. The cropped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

  Figure . we show four synthesis examples. For the first two microtextures (top row in Figure .) one can observe that the generated image does not respect the global statistic of the input creating "stains" of different illuminations.

Figure 5 .

 5 Figure 5.12 -Four crops of different parts of Figure 5.5 at its original resolution. The cropped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

Figure 5 .

 5 Figure 5.13 -Four crops of different parts of Figure 5.6 at its original resolution. The cropped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

  Figure 5.14 -Four crops of different parts of Figure 5.6 at its original resolution. The croped images are of size 500 × 500 pixels. Each one represents a different texture belonging to a single "big texture".

Figure 5 .

 5 Figure 5.15 -Synthesis results of the Portilla and Simoncelli method [50]. It is satisfactory for many small grain textures (first row) but may miss the global structure (second row).

  Figure . and Figure . one can observe these results.

Figure 5 .

 5 Figure 5.16 -Synthesis results of the RPN method [21]. This method works extremely well for a class of "microtextures" including tissues and granular textures with no geometric structures [20].

Figure 5 .

 5 Figure 5.17 -Overview of the Efros and Leung algorithm [14]. Given a texture image (left) a new image (right) is being synthesized a pixel at a time. For a pixel p being synthesized the method finds all neighbourhoods in the left image that match the neighbourhood of p and then chooses randomly one of the neighbourhoods and assigns its central pixel value to p. The image was taken from the work in [14].

Figure 5 .

 5 Figure 5.18 -Synthesis results of the Efros and Freeman method [15]. Works for microtextures with a risk of not respecting the global statistics (see the first experiment, top left). Works for macrotextures with the risk of verbatim copy.

Figure 5 .

 5 Figure 5.20 -Synthesis results of LG [8]. Works for microtextures with a risk of no respect of the global statistics (see the first experiment, top left). Works for macrotextures with a lower risk of verbatim copy than EF for this kind of textures.

Figure 5 Figure 5 Figure 5 Figure 5 Figure 5 Figure 5 Figure 5 Figure 5

 55555555 Figure 5.21 -1) RPN loses line structures. PS homogenizes the input. EF verbatim copy. LG satisfying for this input. 2) RPN loses line structures. PS homogenizes the input. EF satisfying for this input. LG garbage growing. 3) RPN loses line structures. PS satisfying for this input. EF satisfying for this input.LG garbage growing. 4) RPN loses line structures. PS homogenizes the input. EF garbage growing. LG garbage growing. 5) RPN loses line structures. PS homogenizes the input and color artifact. EF verbatim copy and garbage growing. LG verbatim copy.

Figure 5 .

 5 Figure 5.30 -One dimensional anamorphoses functions used to lightly deform the texture images aiming at masking the evident salient structures.

Figure 5 .

 5 Figure 5.35 -Two dimensional anamorphism applied to two different wood textures. The deformed images remain very similar to the input.

Figure 5 .

 5 Figure 5.36 -Results applying the one dimensional random mosaic algorithm to stone and metal textures. From left to right: input, synthesis 1, synthesis 2 and synthesis 3.

Figure 5 .

 5 Figure5.37 -Results applying the one dimensional random mosaic algorithm to a wood texture. From top to bottom: input, synthesis 1, synthesis 2 and synthesis 3.



Figure 5 .

 5 Figure 5.38 -Results applying the one dimensional random mosaic algorithm to a wood texture. From top to bottom: input, synthesis 1, synthesis 2 and synthesis 3.

Figure 5 .

 5 Figure 5.39 -Results applying the one dimensional random mosaic algorithm to a wood texture. From top to bottom: input, synthesis 1, synthesis 2 and synthesis 3.

Figure 5 .

 5 Figure 5.41 -Results applying the two dimensional random mosaic algorithm to stone and metal textures. From left to right: input, synthesis 1, synthesis 2 and synthesis 3.



Figure 5 .

 5 Figure 5.42 -Synthesis results after randomizing the low frequencies of the input sample. From left to right: input, synthesis 1, synthesis 2 and synthesis 3. A strong correlation exists between the low and high frequencies of these input images. Observing these results from far reduces this correlation.

  Figure 5.43 -Texture mixing results. Three different examples are shown. Each example corresponds to three images: inputs I 1 and I 2 and I resulting from mixing I 1 and I 2 .

Figure 5 .

 5 Figure 5.44 -Texture mixing results. Three different examples are shown. Each example corresponds to three images: inputs I 1 and I 2 and I resulting from mixing I 1 and I 2 .

  Figure A.1 -Three different iterations of the synthesis process are shown. At each iteration a patch is being synthesized. This patch is represented by the pink square for the three cases.

  Figure A.3 -Three overlap cases arises in the raster scan order: vertical overlap for the first row, horizontal overlap for the first column, and L-shaped overlap everywhere else.

Algorithm  :

 : Patch Selection  Input: Sample texture I 0 , patch under construction P old , patch size w p , tolerance ε > 0, binary weight Q defining the overlap region  Output: Patch position (m, n) : Compute the squared distance image D (see Equation (A.)) containing the distances between the patch P old and all patches of I 0 . : d min ← min (k,l) D(k, l).

:

  Uniformly draw a patch position (m, n) among the set {(k, l), D(k, l) < (1 + ε)d min }.

PQQ=

  ext (m, n) = P old (m, n) if (m, n) ∈ {0, . . . , w p -1} × {0, . . . , w pj)(P old (i, j) -I 0 (m + i, n + j)) ext (k, l)(P ext (k, l) -I 0 (m + k, n + l)) ext (k, l)P ext (k, l) 2 -2Q ext (k, l)P ext (k, l)I 0 (m + k, n + l) + Q ext (k, l)I 0 (m + k, n + l) 2

Figure

  Figure A.4 -Evolution of the synthesized image with the parallelization. New patches from two subsequent iterations are added simultaneously.



  Input: Sample texture I 0 , N r is the number patch rows, N c is the number of patch columns  Output:Synthesized texture I s : N it ← 2(N r -1) + N c : Initialize I s : for k = 1 to N it -1 do : for i = max(0, ⌊ k-Nc+1 2 ⌋) to min(N r -1, ⌊ k 2 ⌋) do :I s ← synthesize a new patch at position i(w pw o ), (k -2i)(w pw o )

Figure A. 5 -

 5 Figure A.5 -Results representation. From left to right: texture sample, position map, synthesized image and synthesis map. The synthesis map shows for each synthesized patch its position in the texture sample. It allows then to identify exactly the verbatim copy regions (continuous color areas of the map).

  and A. successful results are shown for different type of textures. For each example the patch size w p is adapted taking one of the following values {10, 20, 40, 80}. The remaining parameters are fixed to w o = 0.25 and ε = 0.1. In Figure A. three failure cases are shown. The first one is due to the patch size. When this one is too small, in particular for macro textures, the algorithm fails to recover the details of the different scales as can be seen in the first row examples' in Figure A.. The second failure is the verbatim copy effect. For some texture samples the verbatim copy zones are visually noticeable and unnatural. This is illustrated in the second row examples' in Figure A..

First

  Figure A.6 -Successful results of Efros and Freeman image quilting algorithms. The small images represent the example texture and the big ones the corresponding synthesis result. Each row of examples is followed by a row containing the corresponding color and synthesis maps. For all examples the patch size w p and the ratio r used is indicated. The overlap size is fixed to w o = 0.25 and the tolerance error to ε = 0.1.

  Figure A. for w o = 0.50. On the other hand if w o is low the set of candidate patches is bigger and thus reduces the verbatim copy effect on the cost of decreasing the visual quality of the results as can be seen in the third example in Figure A..

  Figure A.9 -Patch size influence. Each couple of rows shows from left to right the input, the synthesis results for w p = 10, 20, 40, 80 and the corresponding color and synthesis maps. For all examples w o = 0.25 and ε = 0.1.

  Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD  /  Saint-Aubin, France  Title : Exemplar based texture synthesis: models and applications Keywords : exemplar based texture synthesis, Gaussian model, patch, multi scale, optimal transport Abstract : This dissertation contributes to the problem of exemplar based texture synthesis by introducing the use of local Gaussian patch models to generate new texture images.

  L. Raad, A. Desolneux and J.-M. Morel. Can exemplar based texture synthesis models emulate large textures? To be submitted.

	• L. Raad, A. Desolneux and J.-M. Morel. Midway blending as a new technique for patch quilt-ing. To be submitted.
	. Introduction
	Exemplar-based texture synthesis aims at generating new texture images inspired from a texture
	sample. Texture synthesis algorithms can be roughly divided in two categories: the statistics-based
	methods [, , , , ] and the non parametric patch based methods [, , , , , , , ].
	  

• L. Raad, A. Desolneux and J.-M.

Morel. Locally Gaussian Exemplar Based Texture Synthesis.

In IEEE International Conference on Image Processing (ICIP), pages  -, . IEEE. • L. Raad, A. Desolneux and J.-M. Morel. Conditional Gaussian Models for Texture Synthesis. In International Conference on Scale Space and Variational Methods in Computer Vision (SSVM), pages  -, . Springer International Publishing. • L. Raad, A. Desolneux and J.-M. Morel. Multiscale Exemplar Based Texture Synthesis by Locally Gaussian Models. In IberoAmerican Congress on Pattern Recognition (CIARP),pages  -, . Springer Iternational Publishing. • L. Raad, A. Desolneux and J.-M. Morel. Conditional and multiscale locally Gaussian models for texture synthesis. In Journal of Mathematical Imaging and Vision (JMIV), pages  -, . • L. Raad and B. Galerne. Efros and Freeman Image Quilting Algorithm for Texture Synthesis. Submitted to Image Processing On Line (IPOL).

•

Table 2 .

 2 1 -Summary of the notations used in this chapter.

Table 3

 3

.1 -Summary of the notations used in this chapter.

  Simple illustrative example 2. Left image: signals s 1 and s 2 to register. Middle image: cumulative histograms H 1 and H 2 and the midway histogram H. Right image: signals s 1 s 2 and the midway signal s. In this example one can observe that the midway signal is not the local translation of s 1 and s 2 that one would expect. The middle image shows that the resulting midway histogram H (i.e. the harmonic mean of H 1 and H 2 ) has an additional slope with respect to H 1 and H 2 . This means that the midway signal s has an additional intermediary flat zone.

	1.2		2.2
		H1		s1
		H2		s2
	1	Hmid	2	smid
			1.8
	8
			1.6
			1.4
			1.2
			1
			0.8
			0	10	20	30	40	50	60	70	80	90	100
	Figure 4.6 -

  Signal to register s 1 and s 2  Output: Midway signal s mid : n ← length(s 1 ) {equal to length(s 2 )}

	Algorithm : Midway registration of two signals
	}	1), 2/(n -1), . . . , 1}) {see Algorithm
	: return H mid
	

 Input:

Algorithm : Minimum Error Boundary Cut

 Input: Output texture I s , patch P old (from I s ), patch P in (from I 0 ), patch size w p , overlap size w o  Output: Boundary cut γ : Compute the squared difference e(i, j) = (P new (i, j) -P old (i, j)) 2 : switch (overlap type) : case vertical:

Compute the minimal cumulative vertical error E v :

:

E v (w p -1, j) = e(w p -1, j), j ∈ {0, . . . , w o -1}

:

end for :

Determine j * = arg min j E v (0, j)

:

Trace back the path γ starting at γ wp-1 = (0, j * ) using T v : for i = w p -2 to 0,

Compute the minimal cumulative horizontal error E h :

Trace back the path γ starting at γ wp-1 = (i * , 0) using T h : for j = w p -2 to 0,

Compute the minimal cumulative vertical error E v as in the vertical case :

Compute the minimal cumulative horizontal error E h as in the horizontal case :

Trace the vertical part of γ starting at (i * , i * ) using T v : for i = w pi * -2 to 0,

Trace the horizontal part of γ starting at (i * , i * ) using

A.. Blending along the Cut

This is the last step of an iteration. Its goal is to construct the new patch P new by blending P old and P in using the previously computed boundary cut. The boundary cut defines a binary mask M that  equals one on the left and/or top of the cut and zero otherwise. The patch P new can be defined as

We noticed that the Matlab implementation of the authors proposes an optional smoothing of the binary mask presumably to avoid noticeable transitions along the cut. We implemented this option but we did not observe noticeable improvements. Since it is not discussed in the original paper we do not use it for the experimental results of this paper.

A. Implementation

A.. Computing Patch Distances with FFT

In this section we explicit an algorithm to compute the squared distance between a patch P old and all the patches of the input texture I 0 using Fast Fourier Transform (FFT), that is, to compute the squared distance image D of Equation (A.) involved in the patch search algorithm (see Algorithm ). To the best of our knowledge, regarding texture synthesis this acceleration was first discussed by Kwatra et al. [?]. Recall that we consider a patch P old of size w p × w p and that I 0 is of size M 0 × N 0 . Hence the naive summation to compute

p M 0 N 0 operations with w p = 40 pixels as a typical value. An important asset of the proposed FFT-based implementation is that the computational cost is limited to two FFT calls for images of size M 0 × N 0 and is thus independent of the patch width w p .

Let us first recall some notation. A patch of size w p × w p within an image is represented by its top-left corner, hence all admissible patches of I 0 have (m, n)-coordinates in the set {0, . . . , M 0w p } × {0, . . . , N 0w p }.

Discrete Fourier Transform

Convolution Product Denote by V * W the convolution product between V and W , that is,

where the indexes are (mk, nl) are understood modulo (M 0 , N 0 ). With these conventions for the DFT, one has F(V * W ) = M 0 N 0 V Ŵ , where the multiplication between images is the componentwise product. However, recall that the FFTW library computes the Fourier transform