\documentclass{ipol}

\ipolSetTitle{Efros and Freeman Image Quilting Algorithm for Texture Synthesis}
\ipolSetAuthors{Lara Raad\ipolAuthorMark{1},
                Bruno Galerne\ipolAuthorMark{2}}
\ipolSetAffiliations{%
\ipolAuthorMark{1} CMLA, ENS Cachan, France
                   (\texttt{raad@cmla.ens-cachan.fr})\\
\ipolAuthorMark{2} Laboratoire MAP5 (UMR CNRS 8145), Universit\'{e} Paris Descartes, Sorbonne Paris Cit\'{e} (\texttt{bruno.galerne@parisdescartes.fr})	}

\ipolSetID{171}
\ipolSetVolume{6}

\ipolSetSubmissionYear{2016}
\ipolSetSubmissionMonth{02}
\ipolSetSubmissionDay{15}

\ipolSetAcceptedYear{2016}
\ipolSetAcceptedMonth{09}
\ipolSetAcceptedDay{16}

\ipolSetPublicationYear{2016}
\ipolSetPublicationMonth{12}
\ipolSetPublicationDay{00}

\setcounter{page}{1}

\usepackage{hyperref,graphicx,amsmath,amssymb,amsthm}
\usepackage{url}

\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{subfig}
% macros
\newcommand{\R}{\mathbb{R}} % real numbers
\newcommand{\F}{\mathcal{F}} % Fourier transform
\newcommand{\argmin}{\operatorname*{argmin}} % real numbers


% for switch case statement
\newcommand{\SWITCH}[1]{\STATE \textbf{switch} (#1)}
\newcommand{\ENDSWITCH}{\STATE \textbf{end switch}}
\newcommand{\CASE}[1]{\STATE \textbf{case} #1\textbf{:} \begin{ALC@g}}
\newcommand{\ENDCASE}{\end{ALC@g}}
\newcommand{\CASELINE}[1]{\STATE \textbf{case} #1\textbf{:} }
\newcommand{\DEFAULT}{\STATE \textbf{default:} \begin{ALC@g}}
\newcommand{\ENDDEFAULT}{\end{ALC@g}}
\newcommand{\DEFAULTLINE}[1]{\STATE \textbf{default:} }

\newcommand{\Pold}{P_{\text{old}}}
\newcommand{\Pnew}{P_{\text{new}}}
\newcommand{\Pin}{P_{\text{in}}}



%-------------------------------------------------------------------------------
\begin{document}

%-------------------------------------------------------------------------------
\begin{ipolAbstract}
  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 programming problem used for the quilting step as well as implementation details. Additionally we propose a partially parallel version of the quilting technique.
\end{ipolAbstract}

\begin{ipolCode}
The source code and an online demonstration of the algorithm described in this article are accessible at the \href{\ipolLink}{IPOL web page of this article}. 
\end{ipolCode}

%-------------------------------------------------------------------------------
\ipolKeywords{texture synthesis; quilting; patch; linear programming}

%-------------------------------------------------------------------------------
\section{Introduction}


Texture synthesis is a classical image processing problem that finds its applications in virtual reality rendering (video games, animation movies, \dots).
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~\cite{Tartavel_sparsity_spectrum_constraints_2014, Raad_Desolneux_Morel_locally_Gaussian_texture_synthesis_2014}.

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~\cite{Heeger_Bergen_pyramid_based_texture_analysis_synthesis_1995, BVGR_Heeger_Bergen_2014},
several methods are based on statistics of wavelet coefficients or more involved multiscale image representations~\cite{Portilla_Simoncelli_joint_statistics_complex_wavelet_coefficients_2000, Peyre_texture_synthesis_grouplets_2010, Rabin_Peyre_Delon_Bernot_wasserstein_barycenter_texture_mixing_proc_2010}.
Another approach initially proposed by van Wijk~\cite{Wijk_spot_noise_texture_synthesis_1991}, which consists in randomizing the Fourier phase of an image, has been extended to an exemplar-based synthesis method in~\cite{Galerne_Gousseau_Morel_random_phase_textures_2011, Galerne_Gousseau_Morel_micro_texture_synthesis_IPOL_2011}. This method yields excellent results for micro-textures, which are defined as 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~\cite{Efros_Leung_texture_synthesis_non_parametric_sampling_1999, Wei_Levoy_fast_texture_synthesis_2000}.
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~\cite{Aguerrebere_Gousseau_Tartavel_IPOL_Efros_Leung_2013}.
Along with the method of Liang et al.\cite{Liang_et_al_realtime_texture_synthesis_patch_based_2001} developed at the same period,
image quilting~\cite{Efros_Freeman_image_quilting_2001},
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~\ref{sec_algorithm_description} for details).
Let us note that another solution of this ``stitching step''  is proposed in~\cite{Kwatra_et_al_graphcut_textures_2003} 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~\cite{WLKT_soa_examplar_based_texture_2009} for a more complete survey of this category of texture synthesis algorithms. Let us also mention the recent paper~\cite{self_tuning_texture_opt_2015} that discusses and attempts to solve some limitations of these approaches.

The plan of the paper is as follows.
Section~\ref{sec_algorithm_description} describes in detail all the steps of the algorithm and in particular gives a detailed mathematical description of the linear programming problem for the computation of the minimum error boundary cut.
Section~\ref{sec_implementation} gives implementation details and discusses a parallelization of the Efros-Freeman algorithm which, to the best of our knowledge, is a contribution of this paper.
Section~\ref{sec_experiments} presents numerous experiments of our implementation of the Efros-Freeman algorithm.
This experimental section shows that this algorithm usually produces 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~\cite{Efros_Leung_texture_synthesis_non_parametric_sampling_1999, Aguerrebere_Gousseau_Tartavel_IPOL_Efros_Leung_2013}, are also present in the Efros-Freeman results, but only at a larger scale.

\section{Algorithm Description}
\label{sec_algorithm_description}

In~\cite{Efros_Freeman_image_quilting_2001} 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_\text{old}$ of $I_s$ that is only partially defined on a region called overlap region (see Figure~\ref{fig:sequential-alg} and Figure~\ref{fig:overlap_cases}). To do so a patch $P_\text{in}$ of $I_0$ that matches $P_\text{old}$ on the overlap region is randomly selected. An optimal boundary cut between $P_\text{old}$ and $P_\text{in}$ is then computed within the overlap region. This optimal boundary cut is used to construct the new patch $P_\text{new}$ by blending $P_\text{old}$ and $P_\text{in}$ along the cut (see Figure~\ref{fig:quilting_illustration}).

The whole image quilting algorithm is described in Algorithm~\ref{alg:quilting} and the remaining of this section will detail each step of the algorithm, namely the initialization, the patch search procedure to select $P_\text{in}$, the computation of the minimum error boundary and the blending procedure along the boundary.

\begin{figure}[!htb]
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[input]{\includegraphics[width=0.12\linewidth]{figs/input/simoncelli_8.png}}
\quad
\subfloat[iteration 10]{\includegraphics[width=0.25\linewidth]{figs/simoncelli_8_it10.png}}
\quad
\subfloat[iteration 87]{\includegraphics[width=0.25\linewidth]{figs/simoncelli_8_it87.png}}
\quad
\subfloat[iteration 217]{\includegraphics[width=0.25\linewidth]{figs/simoncelli_8_it217.png}}
\caption{Three different iterations of the synthesis process are shown. At each iteration a patch is being synthesized. This patch is represented by the red square for the three cases.}
\label{fig:sequential-alg}
\end{figure}

\begin{figure}[!htb]
\centering
\subfloat[Vertical overlap]{\includegraphics[width=0.3\linewidth]{figs/scan_and_overlap/vertical_overlap.png}}
\quad
\subfloat[Horizontal overlap]{\includegraphics[width=0.3\linewidth]{figs/scan_and_overlap/horizontal_overlap.png}}
\quad
\subfloat[L-shaped overlap]{\includegraphics[width=0.3\linewidth]{figs/scan_and_overlap/Lshaped_overlap.png}}
\caption{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.}
\label{fig:overlap_cases}
\end{figure}

\begin{figure}[!htb]
\centering

\subfloat[]{\includegraphics[width=0.26\linewidth]{figs/patch_1.png}}
\quad
\subfloat[]{\includegraphics[width=0.15\linewidth]{figs/patch_2.png}}
\quad
\subfloat[]{\includegraphics[width=0.15\linewidth]{figs/error_surface.png}}
\quad
\subfloat[]{\includegraphics[width=0.26\linewidth]{figs/quilted_patches.png}}

\caption{Quilting a square patch from the input texture into the synthesized texture. (a) sub-part of the synthesized texture. The red dotted zone shows where the new patch will be quilted. (b) patch to quilt in the red dotted zone shown in (a). (c) error surface between the patch in (b) and the red dotted zone in (a). The red path shows the minimum error boundary. (d) the patch in (b) is quilted along the minimum error boundary in the corresponding zone showed in (a).}
\label{fig:quilting_illustration}
\end{figure}

\begin{algorithm}[ht!]
\caption{Image Quilting}
\label{alg:quilting}
\textbf{Input:} Sample texture $I_{0}$, patch size $w_p$, overlap size $w_o$, tolerance parameter $\varepsilon$, output/input size ratio $r$\\
\textbf{Output:} Synthesized texture $I_{s}$
\begin{algorithmic}[1]
\STATE Initialize $I_s$
\FOR{\textbf{each} patch $P_\text{old}$ in $I_s$}
\STATE Select a compatible patch $P_\text{in}\in I_0$ using the Patch Selection algorithm (see Algorithm~\ref{alg:patch_selection})
\STATE Compute minimum error boundary cut between $P_\text{old}$ and $P_\text{in}$ (see Algorithm~\ref{alg:boundCut})
\STATE Construct the patch $P_\text{new}$ by blending $P_\text{old}$ and $P_\text{in}$ along the boundary cut (see Equation~\eqref{eq:patchCombination})
\STATE Replace $P_\text{old}$ with $P_\text{new}$ within $I_s$
\ENDFOR
\end{algorithmic}
\end{algorithm}

\subsection{Initialization}

The first step of the algorithm is to initialize $I_s$. For that a random patch $P_\text{in}$ of size $w_p\times w_p$ is taken from $I_0$ and placed at the top-left corner of $I_s$.

\subsection{Patch Selection}
\label{sec:PatchSelection}

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_\text{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 for the first row, horizontal overlap for the first column, and L-shaped overlap everywhere else (see Figure~\ref{fig:overlap_cases}).

To select a patch $P_\text{in}$ of the input image $I_0$ one computes the square distance between the overlap region of the patch $P_\text{old}$ of $I_s$ and the corresponding regions of all the patches of $I_0$. The minimal distance $d_\text{min}$ is determined and $P_\text{in}$ is randomly picked among all patches whose distance to $P_\text{old}$ is lower than $(1+\varepsilon)d_\text{min}$ where $\varepsilon$ 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)\in\{ 0,\dots, M_0-w_p \} \times \{ 0,\dots, N_0-w_p \}$.
The squared distance image $D$ contains at each position $(m,n)$ the distance between $P_\text{old}$
and the patch from $I_0$ whose 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)\in\{ 0,\dots, M_0-w_p \} \times \{ 0,\dots, N_0-w_P \}$, one has
\begin{equation}
  D(m,n) = \sum_{i,j=0}^{w_p-1} Q(i,j)(P_\text{old}(i,j) - I_0(m+i,n+j))^2.
\label{eq:distance}
\end{equation}
The computation of the squared distance image $D$ is discussed in Section~\ref{sec_implementation_distances} and the whole patch selection procedure is summarized in Algorithm~\ref{alg:patch_selection}.

\begin{algorithm}[ht!]
\caption{Patch Selection}
\label{alg:patch_selection}
\textbf{Input:} Sample texture $I_{0}$, patch under construction $P_\text{old}$, patch size $w_p$, tolerance $\varepsilon>0$, binary weight $Q$ defining the overlap region\\
\textbf{Output:} Patch position $(m,n)$\\
\begin{algorithmic}[1]
\STATE Compute the squared distance image $D$ (see Equation~\eqref{eq:distance}) containing the distances between the patch $P_\text{old}$ and all patches of $I_0$.
\STATE $\displaystyle d_\text{min} \leftarrow \min_{(k,l)} D(k,l)$.
\STATE Uniformly draw a patch position $(m,n)$ among the set $\{(k,l),~D(k,l) < (1+\varepsilon)d_\text{min}\}$.
\end{algorithmic}

\end{algorithm}

\subsection{Minimum Error Boundary Cut}
\label{subsec:minErrorCut}

This step of the algorithm is the main contribution of the Efros-Freeman algorithm~\cite{Efros_Freeman_image_quilting_2001}.

The Patch Search step (see Section~\ref{sec:PatchSelection}) gives a patch $P_\text{in}$ of $I_0$ having coordinates $(m,n)$
that is similar to the partially defined current patch $P_\text{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~\ref{fig:overlap_cases}.
To get the final patch $P$ one must combine the patches $P_\text{old}$ and $P_\text{in}$.
Denoting $Q$ the binary weight for the overlap regions as in the previous section, then, for any binary image $M$ such that $0\leq M(i,j)\leq Q(i,j),~(i,j)\in\{1,\dots,w_p\}^2$, $P$ can be defined as the combination
\begin{equation*}
\Pnew = M P_{\text{old}} + \left(1-M\right)P_{\text{in}}.
\end{equation*}
The main idea of Efros and Freeman~\cite{Efros_Freeman_image_quilting_2001} is looking for a binary shape $M$ where the transition between $P_{\text{old}}$ and $P_{\text{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 detail the computation of the minimum error boundary cut for the three overlap cases, starting with the vertical and horizontal cases.

\subsubsection{Connected Paths and Boundary Error}

In the following, we call a \textit{path} of length $\ell\geq 1$ in $\{0,\dots,w_p-1\}^2$ any ordered sequence $\gamma = (\gamma_0,\gamma_1,\dots, \gamma_{\ell-1})$ of 8-connected pixels $\gamma_k\in\{0,\dots,w_p-1\}^2$, that is,
for all $k\in\{0,\dots,\ell-2\}$, $\max(|\gamma_{k+1}^1 - \gamma_k^1|,|\gamma_{k+1}^2 - \gamma_k^2|) = 1$
(for convenience, within this section the first and second coordinates of a pixel $\gamma_k$ are denoted by $\gamma_k^1$ and $\gamma_k^2$).

Let us recall that we use matrix coordinates for the pixels.
Hence a path $\gamma$ of length $\ell$ is said to be
\textit{vertical} if for all $k\in\{0,\dots,\ell-2\}$,
$\gamma_{k+1}^1 - \gamma_k^1 = -1$ (vertical paths are oriented from bottom to top),
and \textit{horizontal} if for all $k\in\{0,\dots,\ell-2\}$,
$\gamma_{k+1}^2 - \gamma_k^2 = -1$ (horizontal paths are oriented from right to left).

Denoting by $e(i,j) = (P_{\text{new}}(i,j) - P_{\text{old}}(i,j))^2$ the squared difference between the two patches $P_{\text{old}}$ and $P_{\text{new}}$,
the boundary error of a path $\gamma$ of length $\ell$ is defined by
$$
\mathcal{E}(\gamma) = \sum_{k=0}^{\ell-1} e(\gamma_k).
$$


\subsubsection{Vertical Boundary Cuts}

In this section we explain how the optimal boundary between $P_{\text{old}}$ and $P_{\text{new}}$ is defined and computed in the case of vertical overlap, that is when the overlap region is the rectangle $\{0,\dots,w_p-1\}\times\{0,\dots,w_o-1\}$.

The optimal boundary is defined as the vertical path $\gamma$ that minimizes the boundary error $\mathcal{E}(\gamma)$ while joining both ends of the overlap regions.
More precisely, define the admissible vertical paths as
$$
\Gamma_v
= \{ \gamma=(\gamma_0,\dots,\gamma_{w_p-1}),~\forall k,~\gamma_k^1 = w_p-1-k ~\text{and}~\gamma_k^2\in\{0,\dots,w_o-1\}\},
$$
(one notices that paths of $\gamma\in\Gamma_v$ are indeed vertical since $\gamma_{k+1}^1 - \gamma_k^1 = -1$).
Then an optimal boundary is defined as any solution of the optimization problem
\begin{equation}
\min_{\gamma\in\Gamma_v} \mathcal{E}(\gamma).
\label{eq:problem_vertical_cut}
\end{equation}

Problem~\eqref{eq:problem_vertical_cut} is solved using dynamic programming.
This is possible because Problem~\eqref{eq:problem_vertical_cut} verifies the principle of optimality:
if ${(\gamma_0, \dots, \gamma_{w_p-1})}$ is an optimal solution of Problem~\eqref{eq:problem_vertical_cut}, 
then, for all $0\leq r\leq w_p-1$, the subpath $(\gamma_0, \dots, \gamma_{r-1})$ is an optimal vertical path to reach the $r$-th point $\gamma_{r-1}$ starting from the bottom of the overlap region.

Let us now discuss in detail the dynamic programming method to solve Problem~\eqref{eq:problem_vertical_cut}.
For all points $(i,j)\in\{0,\dots,w_p-1\}\times\{0,\dots,w_o-1\}$, denote by
$\Gamma_v^{(i,j)}$
the set of vertical paths $\gamma$ that start at the bottom side of the overlap region and end at the point $(i,j)$, that is
$$
\Gamma_v^{(i,j)}
= \{ \gamma=(\gamma_0,\dots,\gamma_{w_p-1-i}),~\forall k,~\gamma_k^1 = w_p-1-k,~\gamma_k^2\in\{0,\dots,w_o-1\},~\text{and}~\gamma_{w_p-1-i}=(i,j) \}.
$$
Now,  for all points $(i,j)\in\{0,\dots,w_p-1\}\times\{0,\dots,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_{\gamma\in\Gamma_v^{(i,j)}}  \mathcal{E}(\gamma).
$$
Then, since
$$
\Gamma_v = \bigcup_{j\in\{0,\dots,w_o-1\}} \Gamma_v^{(0,j)},
$$
one has
$$
\min_{\gamma\in\Gamma_v} \mathcal{E}(\gamma)
= \min_{j\in\{0,\dots,w_o-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\in\{0,\dots,w_p-2\}$, one has
\begin{equation}
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) ),
\label{eq:recursive_relation_Evij}
\end{equation}
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~\eqref{eq:problem_vertical_cut} is to compute the costs $E_v(i,j)$ line by line from bottom to top using Equation~\eqref{eq:recursive_relation_Evij} recursively, and then search for the minimal value of the first line $j^* = \argmin_j E_v(0,j)$.
The full optimal path $\gamma$ can then be traced back starting from this coordinate $(0,j^*) = (\gamma^1_{w_p-1}, \gamma^2_{w_p-1})$.
More precisely, if $(i,j)$ are the coordinates of the point $\gamma_{w_p-1-i}$ of the optimal path $\gamma$, then its preceding point $\gamma_{w_p-1-(i+1)}$ is given by
$$
\gamma_{w_p-1-(i+1)} = \argmin( 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) = \argmin( 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~\eqref{eq:recursive_relation_Evij} in order to trace back the path without additional computation since for $i=w_p - 2$ to $0$,  $\gamma_{i} = T_v(\gamma_{i+1})$.


\subsubsection{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
$$
\Gamma_h
= \{ \gamma=(\gamma_0,\dots,\gamma_{w_p-1}),~\forall k,~\gamma_k^2 = w_p-1-k ~\text{and}~\gamma_k^1\in\{0,\dots,w_o-1\}\},
$$
and for all points $(i,j)\in\{0,\dots,w_o-1\}\times\{0,\dots,w_p-1\}$,
one defines
$$
\Gamma_h^{(i,j)}
= \{ \gamma=(\gamma_0,\dots,\gamma_{w_p-1-j}),~\forall k,~\gamma_k^2 = w_p-1-k,~\gamma_k^1\in\{0,\dots,w_o-1\},~\text{and}~\gamma_{w_p-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_{\gamma\in\Gamma_h^{(i,j)}}  \mathcal{E}(\gamma),
$$
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\in\{0,\dots,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)\in\{0,\dots,w_o-1\}\times\{0,\dots,w_p-1\}$.
Since
$$
\min_{\gamma\in\Gamma_h}  \mathcal{E}(\gamma)
= \min_{i\in\{0,\dots,w_o-1\}} E_h(i,0)
$$
the end point of the optimal horizontal path is $(i^*,0)$ where $i^* = \argmin_i E_h(i,0)$ and
the path can be traced back thanks to the matrix
$T_h(i,j) = \argmin( E_h(i-1,j+1), E_h(i,j+1), E_h(i+1,j+1) )$.

\subsubsection{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~\cite{Efros_Freeman_image_quilting_2001} 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\in\{0,\dots,w_o-1\}$)
and end at the right side of the L-shaped (that is, the points $(i,j)$ with $i\in\{0,\dots,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\in\{0,\dots,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\in\{0,\dots,w_o-1\}$,
one defines
$$
\Gamma_L^i
= 
\{ \gamma = (\gamma_0,\dots, \gamma_{2(w_p - i)-1}),
~(\gamma_0,\dots, \gamma_{w_p - i-1})\in\Gamma_v^{(i,i)}
~\text{and}
~(\gamma_{2(w_p - i)-1},\dots, \gamma_{w_p - i-1})\in\Gamma_h^{(i,i)}
\},
$$
and for all $\gamma\in \Gamma_L^i$ we denote by $\gamma^v$ and $\gamma^h$ its associated vertical and horizontal paths of $\Gamma_v^{(i,i)}$ and $\Gamma_h^{(i,i)}$ respectively.
The set of the admissible L-shaped boundaries $\Gamma_L$ is then defined as the disjoint union of the sets $\Gamma^i_L$, that is,
$$
\Gamma_L = \bigcup_{i=0}^{w_0-1} \Gamma^i_L.
$$
As before, we search for an optimal L-shaped boundary cut $\gamma\in \Gamma_L$ having minimal boundary error
$$
\mathcal{E}(\gamma) = \sum_{k=0}^{\ell(\gamma)-1} e(\gamma_k),
$$
where $\ell(\gamma)$ is the length of the path $\gamma$ (which is equal to $2(w_p - i)-1$ if $\gamma\in \Gamma^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 $\gamma\in\Gamma_L^i$ one has
$$
\mathcal{E}(\gamma) = \mathcal{E}(\gamma^v) + \mathcal{E}(\gamma^h) - e(i,i).
$$
Hence one has
$$
\begin{aligned}
\min_{\gamma \in \Gamma_L} \mathcal{E}(\gamma)
&= \min_{i\in\{0,\dots,w_o-1\}} \min_{\gamma \in \Gamma^i_L} \mathcal{E}(\gamma)
\\
&= \min_{i\in\{0,\dots,w_o-1\}} \min_{\gamma \in \Gamma^i_L} \mathcal{E}(\gamma^h) + \mathcal{E}(\gamma^v) - e(i,i)
\\
&=
\min_{i\in\{0,\dots,w_o-1\}}
\left[
\left( \min_{\gamma^v\in\Gamma_v^{(i,i)}} \mathcal{E}(\gamma^v)\right)
+
\left( \min_{\gamma^h\in\Gamma_h^{(i,i)}} \mathcal{E}(\gamma^h)\right)
-
e(i,i)
\right]
\\
& = \min_{i\in\{0,\dots,w_o-1\}} E_v(i,i) + E_h(i,i) - e(i,i),
\end{aligned}
$$
where $E_v$ and $E_h$ are the vertical and horizontal cumulative errors defined in the previous sections.
Hence,
\begin{equation*}
\min_{\gamma \in \Gamma_L} \mathcal{E}(\gamma)  = \min_{i\in\{0,\dots,w_o-1\}} E_v(i,i) + E_h(i,i) - e(i,i).
\end{equation*}
The above equation enables us to determine the optimal position $(i^*,i^*)$ on the diagonal of the optimal path $\gamma$ once the matrices $E_v$ and $E_h$ have been computed recursively.
We can then trace back the vertical part $\gamma^v$ and the horizontal part $\gamma^h$ of the optimal path $\gamma$ using the matrices $T_v$ and $T_h$.
Let us remark that the ``overall minimum'' evoked in~\cite{Efros_Freeman_image_quilting_2001} 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)$.

\begin{algorithm}[ht!]
\caption{Minimum Error Boundary Cut}
\label{alg:boundCut}
\textbf{Input:} Output texture $I_s$, patch $\Pold$ (from $I_s$), patch $\Pin$ (from $I_0$), patch size $w_p$, overlap size $w_o$ \\
\textbf{Output:} Boundary cut $\gamma$\\

\begin{algorithmic}[1]
  \STATE Compute the squared difference $e(i,j) = (\Pin(i,j) - \Pold(i,j))^2$
\SWITCH {overlap type}
\CASE {vertical}
%\IF{vertical overlap}
\STATE Compute the minimal cumulative vertical error $E_v$:
\STATE $E_v(w_p-1,j) = e(w_p-1,j), \quad j\in \{0,\dots, w_o-1 \}$
\FOR{$i= w_p -2$ to $0$}
\STATE $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)), \quad j\in \{0,\dots, w_o-1 \}$
\STATE $T_v(i,j) = \argmin (E_v(i+1,j-1), E_v(i+1,j), E_v(i+1,j+1)), \quad j\in \{0,\dots, w_o-1 \}$
\ENDFOR
\STATE Determine $ j^* = \argmin_j E_v(0,j)$
\STATE Trace back the path $\gamma$ starting at $\gamma_{w_p - 1} = (0,j^{*})$ using $T_v$: for $i=w_p - 2$ to $0$,  $\gamma_{i} = T_v(\gamma_{i+1})$
\ENDCASE
%\ELSIF{horizontal overlap}
\CASE {horizontal}
\STATE Compute the minimal cumulative horizontal error $E_h$:
\STATE $E_h(i,w_p-1) = e(i,w_p-1), \quad i\in \{0,\dots, w_o-1 \}$
\FOR{$j= w_p -2$ to $0$}
\STATE $E_h(i,j) = e(i,j) +\min(E_h(i-1,j+1), E_h(i,j+1), E_h(i+1,j+1)), \quad i\in \{0,\dots, w_o-1 \}$
\STATE $T_h(i,j) = \argmin (E_h(i-1,j+1), E_h(i,j+1), E_h(i+1,j+1)), \quad i\in \{0,\dots, w_o-1 \}$
\ENDFOR
\STATE Determine $ i^* = \argmin_{i} E_h(i,0)$
\STATE Trace back the path $\gamma$ starting at $\gamma_{w_p - 1} = (i^{*},0)$ using $T_h$: for $j=w_p - 2$ to $0$, $\gamma_{j} = T_h(\gamma_{j+1})$
\ENDCASE
\CASE {L-shaped}
%\ELSIF{L-shaped overlap}
\STATE Compute the minimal cumulative vertical error $E_v$ as in the vertical case
\STATE Compute the minimal cumulative horizontal error $E_h$ as in the horizontal case
\STATE Determine $ i^* = \argmin_{i} (E_v(i,i) +E_h(i,i)-e(i,i))$ (with $i \in \{0, \dots, w_o-1\}$)
\STATE Trace the vertical part of $\gamma$ starting at $(i^{*},i^*)$ using $T_v$: for $i=w_p-i^*-2$ to $0$, $\gamma_i = T_v(\gamma_{i+1})$
\STATE Trace the horizontal part of $\gamma$ starting at $(i^{*},i^*)$ using $T_h$: for $j=w_p-i^*$ to $2(w_p-i^*)-1$, $\gamma_j = T_h(\gamma_{j-1})$
\ENDCASE
%\ENDIF
\ENDSWITCH
\end{algorithmic}

\end{algorithm}

\subsection{Blending along the Cut}

This is the last step of an iteration. Its goal is to construct the new patch $P_\text{new}$ by blending $P_\text{old}$ and $P_\text{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 $\Pnew$ can be defined as

\begin{equation}
\Pnew = M P_{\text{old}} + \left(1-M\right)P_{\text{in}}.
\label{eq:patchCombination}
\end{equation}

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 rough way of doing this is considering that all pixels to the left and/or top of the cut correspond to the already synthesized part and all the rest to the new patch. Considering the binary image $M$ defined in section~\ref{subsec:minErrorCut}, the patch under construction $P$ is then obtained as the combination of the old and new information as follows as in~\eqref{eq:patchCombination}.
%
%However there are textures for which the transitions along the cut are slightly noticeable. The authors propose to avoid this artifact by smoothing the transition by applying a smoothing kernel to the binary mask $M$. We used a Gaussian kernel $G_$
%
%To avoid this the authors add an extra step to smooth transitions that consists in Instead of using directly the binary image $M$ they apply to it a smoothing filter, for example a Gaussian kernel $G_{\sigma}$ and obtain like this a smoother image $\tilde{M}$ (no longer binary)
%$$\tilde{M} = M\ast G_{\sigma}$$
%and the new patch $P$ is obtained as in (\ref{eq:blendpatch}). This is illustrated in Figure.\ref{fig:smoothing_illustration}.
%\begin{equation}
%\label{eq:blendpatch}
%P = \tilde{M}P_{\text{old}} + \left(1-\tilde{M}\right)P_{\text{new}}.
%\end{equation}
%
%\begin{figure}[]
%\centering
%
%\subfloat[$P_{\text{old}}$]{\includegraphics[width=0.18\linewidth]{figs/p_old.png}}
%\quad
%\subfloat[$\tilde M$]{\includegraphics[width=0.18\linewidth]{figs/m_tilde.png}}
%\quad
%\subfloat[$P_{\text{new}}$]{\includegraphics[width=0.18\linewidth]{figs/p_new.png}}
%\quad
%\subfloat[$1- \tilde M$]{\includegraphics[width=0.18\linewidth]{figs/m_tilde_inv.png}}
%\quad
%\subfloat[$P$]{\includegraphics[width=0.18\linewidth]{figs/quilted_smooth.png}}
%
%\caption{(a) the patch under construction. (b) the smoothed image $\tilde M$ that determines which part of the old patch is kept. (c) the patch chosen in $I_0$ to blend into $I_s$. (d) the smoothed image $1-\tilde M$ that determines which part of the new patch is kept. (e) the old and new patch blended according to the smoothed images.}
%\label{fig:smoothing_illustration}
%\end{figure}

\section{Implementation}
\label{sec_implementation}

\subsection{Computing Patch Distances with FFT}
\label{sec_implementation_distances}

In this section we explicit an algorithm to compute the squared distance between a patch $\Pold$ and all the patches of the input texture $I_0$ using the Fast Fourier Transform (FFT), that is, to compute the squared distance image $D$ of Equation~\eqref{eq:distance} involved in the patch search algorithm (see Algorithm~\ref{alg:patch_selection}).
To the best of our knowledge, regarding texture synthesis this acceleration was first discussed by Kwatra et al.~\cite{Kwatra_et_al_graphcut_textures_2003}.
Recall that we consider a patch $\Pold$ of size $w_p\times w_p$ and that $I_0$ is of size $M_0\times N_0$.
Hence the naive summation to compute
$$
D(m,n) = \sum_{i,j=0}^{w_p-1} Q(i,j)(\Pold(i,j) - I_0(m+i,n+j))^2,
$$
for all $(m,n)\in \{ 0,\dots, M_0-w_p \} \times \{ 0,\dots, N_0-w_p \}$ requires around $w_p^2 M_0N_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\times N_0$ and is thus independent of the patch width $w_p$.


Let us first recall some notation.
A patch of size $w_p\times 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,\dots, M_0-w_p \} \times \{ 0,\dots, N_0-w_p \}$.

\paragraph{Discrete Fourier Transform.}
Given any image $V\in \R^{M_0\times N_0}$ we denote by $\F(V) = \hat{V}$ the discrete Fourier transform of $V$ and
$\F^{-1}(V) = \check{V}$ the inverse discrete Fourier transform of $V$, defined by
$$
\hat{V}(k,l) = \frac{1}{M_0N_0} \sum_{m=0}^{M_0-1} \sum_{n=0}^{N_0-1} V(m,n) e^{-2i\pi\left(\frac{k m}{M_0} + \frac{n l}{N_0}\right)}
\quad
\text{and}
\quad
\check{V}(k,l) = \sum_{m=0}^{M_0-1} \sum_{n=0}^{N_0-1} V(m,n) e^{2i\pi\left(\frac{k m}{M_0} + \frac{n l}{N_0}\right)}.
$$

\paragraph{Convolution Product.}
Denote by $V\ast W$ the convolution product between $V$ and $W$, that is,
$$
V\ast W(m,n) = \sum_{k=0}^{M_0-1} \sum_{l=0}^{N_0-1} V(k,l) W(m-k, n-l),
$$
where the indexes $(m-k, n-l)$ are understood modulo $(M_0,N_0)$.
With these conventions for the DFT, one has
$\F(V\ast W) = M_0 N_0 \hat{V} \hat{W}$,
where the multiplication between images is the component-wise product.
However, recall that the FFTW library computes the Fourier transform without normalization.
Hence the operators that are computed are respectively $M_0N_0 \F$ for the forward transform and $\F^{-1}$ for the backward transform.
Hence we have
$$
V\ast W = \F^{-1}(\F(V\ast W)) = \frac{1}{M_0N_0} \F^{-1}( M_0N_0 \F(V) M_0N_0 \F(W)).
$$
This means that after multiplying the two FFTW forward transforms and performing the backward inverse transform, one has to normalize dividing by the size of the images.

\paragraph{Cross-correlation.}

Let us denote by $\Gamma(V,W)$ the cross-correlation between two images,
$$
\Gamma(V,W)(m,n) = \sum_{k=0}^{M_0-1}\sum_{l=0}^{N_0-1} V(k,l) W(m+k,n+l).
$$
Note also $\tilde{V}$ the symmetric of $V$ with respect to the origin, that is, $\tilde{V}(m,n) = V(-m,-n)$.
Note that one has
$\Gamma(W,V)(m,n) = \Gamma(V,W)(-m,-n)$,
that is, $\Gamma(W,V) = \widetilde{\Gamma(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,
$$
\Gamma(V,W)(m,n) = \tilde{V}\ast W(m,n) = V\ast \tilde{W} (-m,-n).
$$
Hence the computational cost for a cross-correlation image is the same that for a convolution product, that is three FFT calls.

\smallskip

We need to compute the squared distance between a patch $\Pold$ 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)\in\{ 0,\dots, M_0-w_p \} \times \{ 0,\dots, N_0-w_p \}$, we want to compute
$$
D(m,n) = \sum_{i,j=0}^{w_p-1} Q(i,j)(\Pold(i,j) - I_0(m+i,n+j))^2.
$$
Let us denote by $P_\textrm{ext}$ and $Q_\textrm{ext}$ the extensions of $\Pold$ and $Q$ into images of size $M_0\times N_0$ by filling the domain with $0$-valued pixels, that is
$$
P_\textrm{ext}(m,n)
=
\begin{cases}
\Pold(m,n) & \text{if $(m,n)\in \{ 0,\dots, w_p-1 \} \times \{ 0,\dots, w_p-1 \}$},
\\
0 & \text{otherwise},
\end{cases}
$$
and similarly for $Q_\textrm{ext}$.
Then one has
$$
\begin{aligned}
  D(m,n)  &  = \sum_{i,j=0}^{w_p-1} Q(i,j)(\Pold(i,j) - I_0(m+i,n+j))^2
\\
& = \sum_{k=0}^{M_0-1}\sum_{l=0}^{N_0-1} Q_\textrm{ext}(k,l)(P_\textrm{ext}(k,l) - I_0(m+k,n+l))^2
\\
& = \sum_{k=0}^{M_0-1}\sum_{l=0}^{N_0-1} \left(Q_\textrm{ext}(k,l) P_\textrm{ext}(k,l)^2 - 2 Q_\textrm{ext}(k,l)P_\textrm{ext}(k,l) I_0(m+k,n+l) +  Q_\textrm{ext}(k,l)I_0(m+k,n+l)^2\right)
\\
& = \sum_{i,j=0}^{w_p-1} Q(i,j) \Pold(i,j)^2 - 2 \Gamma(Q_\textrm{ext}P_\textrm{ext}, I_0)(m,n) + \Gamma(Q_\textrm{ext}, I_0^2)(m,n),
\end{aligned}
$$
where the multiplications $Q_\textrm{ext}P_\textrm{ext}$ and $I_0^2$ are component-wise.

Within our procedure, we will compare a patch $\Pold$ 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 $\Pold$ are set to $0$ (by initializing $I_s$ to $0$), one has
$$
\forall (i,j)\in\{0,w_p-1\}^2,\quad Q(i,j) \Pold(i,j) = \Pold(i,j).
$$
Consequently,
$$
\sum_{i,j=0}^{w_p-1} Q(i,j) \Pold(i,j)^2 = \sum_{i,j=0}^{w_p-1} \Pold(i,j)^2 = \|\Pold\|_2^2,
$$
and $Q_\textrm{ext}P_\textrm{ext} = P_\textrm{ext}$, and thus $\Gamma(Q_\textrm{ext}P_\textrm{ext}, I_0) = \Gamma(P_\textrm{ext}, I_0)$.
Hence, the binary mask $Q$ is only influential when computing $\Gamma(Q_\textrm{ext}, I_0^2)$, a computation that is done just once at the beginning of the procedure since it does not depend on $\Pold$.

In the end, the image $D$ is simply given by
\begin{equation}
D = \|\Pold \|_2^2 - 2 \Gamma(P_\textrm{ext}, I_0) + \Gamma(Q_\textrm{ext}, I_0^2),
\label{eq:distFFT}
\end{equation}
where the last image $\Gamma(Q_\textrm{ext}, I_0^2)$ 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 cross-correlation
$\Gamma(P_\textrm{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)\in\{0,\dots,M_0-1\}\times\{0,\dots,N_0-1\}$, but the squared distances for
$(m,n)$-coordinates outside $\{ 0,\dots, M_0-w_p \} \times \{ 0,\dots, N_0-w_p \}$ correspond to patches of size $w_p\times w_p$ that are defined by periodic boundary conditions and thus those patches must be discarded when searching for the minimal distance.

\paragraph{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~\eqref{eq:distFFT} 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~\ref{alg:patch_selection}) since then none of patches satisfy the condition $D(m,n) < (1+\varepsilon) 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$).

\subsection{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 to 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 Figure~\ref{fig:parallelization}.

Let us discuss the acceleration induced by the parallelization. Suppose we synthesize an image $I_s$ made of $N_r\times 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_{\text{it}}$ is $N_{\text{it}} = N_r N_c$.  When parallelizing the number of iterations $N_{\text{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,\left\lfloor \frac{N_c+1}{2} \right\rfloor)}$, where $\left\lfloor x\right\rfloor = \max ( n\in \mathbb{N}, n\leq x)$, since within an iteration there is at most one new patch per row and one new patch per pair of columns.

The parallel version of the quilting algorithm is summarized in Algorithm~\ref{alg:parallelization}. Let us notice that for this parallelization to be valid the overlap size must satisfy $w_o \leq w_p/2$ which is always the case in practice.

\begin{figure}[!htb]
\captionsetup[subfigure]{labelformat=empty}
\centering
\subfloat[iteration $1$]{\includegraphics[width=0.17\linewidth]{figs/T_1-1}}
\hspace{1em}
\subfloat[iteration $2$]{\includegraphics[width=0.17\linewidth]{figs/T_1-2}}
\hspace{1em}
\subfloat[iteration $3$]{\includegraphics[width=0.17\linewidth]{figs/T_1-3}}
\hspace{1em}
\subfloat[iteration $4$]{\includegraphics[width=0.17\linewidth]{figs/T_1-4}}
\hspace{1em}
\subfloat[iteration $5$]{\includegraphics[width=0.17\linewidth]{figs/T_1-5}}

\subfloat[iteration $6$]{\includegraphics[width=0.17\linewidth]{figs/T_1-6}}
\hspace{1em}
\subfloat[iteration $7$]{\includegraphics[width=0.17\linewidth]{figs/T_1-7}}
\hspace{1em}
\subfloat[iteration $8$]{\includegraphics[width=0.17\linewidth]{figs/T_1-8}}
\hspace{1em}
\subfloat[iteration $9$]{\includegraphics[width=0.17\linewidth]{figs/T_1-9}}
\hspace{1em}
%\subfloat[Iteration $10$]{\includegraphics[width=0.17\linewidth]{figs/T_1-10}}
\subfloat[iteration $N_{\text{it}}$]{\includegraphics[width=0.17\linewidth]{figs/T_1-FULL}}

\caption{Evolution of the synthesized image with the parallelization. New patches from two subsequent iterations are added simultaneously.}
\label{fig:parallelization}
\end{figure}

\begin{algorithm}[!htb]
\caption{Parallelization}
\label{alg:parallelization}
\textbf{Input:} Sample texture $I_0$, $N_r$ is the number patch rows, $N_c$ is the number of patch columns\\
\textbf{Output:} Synthesized texture $I_s$\\

\begin{algorithmic}[1]
\STATE $N_{\text{it}}\leftarrow 2(N_r-1) + N_c$
\STATE Initialize $I_s$
\FOR{ $k=1$ \TO $N_{\text{it}}-1$ } 
	\FOR{ $i=\max(0,\lfloor \frac{k-N_c+1}{2}\rfloor)$ \TO $\min(N_r-1,\lfloor \frac{k}{2}\rfloor)$ }
	 	\STATE $I_s \leftarrow $ synthesize a new patch at position $\big(i(w_p-w_o),(k-2i)(w_p-w_o)\big)$
        \rlap{\smash{$\left.\begin{array}{@{}c@{}}\\{}\\{}\end{array}\color{red}\right\}%
          \color{red}\begin{tabular}{l}This loop is\\run in parallel\end{tabular}$}}
	\ENDFOR
\ENDFOR
\end{algorithmic}
\end{algorithm}

\section{Experiments}
\label{sec_experiments}

In this section texture synthesis results are shown using the algorithm described previously~\cite{Efros_Freeman_image_quilting_2001}. Two aspects of this method are emphasized. On the one hand the visual quality of the 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~\ref{fig:result-representation}: 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 shows for each synthesized patch its initial position in the texture sample. It allows to identify exactly the verbatim copy regions which correspond to the continuous color areas of the synthesis map. 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$.

\begin{figure}[!htb]
\captionsetup[subfigure]{labelformat=empty}
\centering

\subfloat{\includegraphics[width=0.15\linewidth]{figs/results_representation/ground1011.png}}
\quad
\subfloat{\includegraphics[width=0.15\linewidth]{figs/results_representation/ground1011_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.3\linewidth]{figs/results_representation/ground1011_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=0.3\linewidth]{figs/results_representation/ground1011_w40_r2_map_out.png}}

\caption{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).}
\label{fig:result-representation}
\end{figure}

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 $\varepsilon$. The overlap size is expressed as the proportion with respect to the patch size. That is $w_o = 0.25$ implicitly means $w_o = 0.25w_p$.

In Figures~\ref{fig:results-success1} and~\ref{fig:results-success2} successful results are shown for different types of textures. For each example the patch size $w_p$ is adapted taking one of the following values $\lbrace 10, 20, 40, 80\rbrace$. The remaining parameters are fixed to $w_o = 0.25$ and $\varepsilon = 0.1$. One can notice that for some texture samples, for example the brick wall in Figure~\ref{fig:results-success1} and the texture samples of the last row in Figure~\ref{fig:results-success2}, even though the results are visually satisfying they are made of large parts of the input sample.

In Figure~\ref{fig:results-failures} three failure cases are shown. The first one is due to the patch size. When this 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 of Figure~\ref{fig:results-failures}. 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 of Figure~\ref{fig:results-failures}. The last failure case is the ``growing garbage'' drawback, which occurs 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~\ref{fig:results-failures} 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.
%
\newlength{\widthfiguresuccessA}
\setlength{\widthfiguresuccessA}{0.096\linewidth}

\begin{figure}[!htbp]
\captionsetup[subfigure]{labelformat=empty}
\centering

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/ground1003.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/ground1003_out_w20_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/IMG_9955.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/IMG_9955_out_w20_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/IMG_9946.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/IMG_9946_out_w40_r2.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/ground1003_map_in.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/ground1003_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/IMG_9955_map_in.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/IMG_9955_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/IMG_9946_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/IMG_9946_w40_r2_map_out.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/IMG_9973.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/IMG_9973_out_w20_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/fabric1003.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/fabric1003_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/grid1004.png}}
\quad
\subfloat[$w_p=80, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/grid1004_out_w80_r2.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/IMG_9973_map_in.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/IMG_9973_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/fabric1003_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/fabric1003_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/grid1004_map_in.png}}
\quad
\subfloat[$w_p=80, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/grid1004_w80_r2_map_out.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/tissus4.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/tissus4_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_11.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_11_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_103.png}}
\quad
\subfloat[$w_p=80, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_103_out_w80_r2.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/tissus4_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/tissus4_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_11_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_11_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_103_map_in.png}}
\quad
\subfloat[$w_p=80, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_103_w80_r2_map_out.png}}

\caption{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 the 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 $\varepsilon = 0.1$.}
\label{fig:results-success1}
\end{figure}

\begin{figure}[!htbp]
\captionsetup[subfigure]{labelformat=empty}
\centering

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/ground1013.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/ground1013_out_w20_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_73.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_73_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_29.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_29_out_w40_r2.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/ground1013_map_in.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/ground1013_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_73_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_73_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/simoncelli_29_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/simoncelli_29_w40_r2_map_out.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/visTex_2.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/visTex_2_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/tissus10.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/tissus10_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/visTex_15.png}}
\quad
\subfloat[$w_p=80, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/visTex_15_out_w80_r2.png}}

\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/visTex_2_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/visTex_2_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/tissus10_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/tissus10_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfiguresuccessA]{figs/input/visTex_15_map_in.png}}
\quad
\subfloat[$w_p=80, r=2$]{\includegraphics[width=2\widthfiguresuccessA]{figs/results_25/output/visTex_15_w80_r2_map_out.png}}

\caption{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 the 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 $\varepsilon = 0.1$.}
\label{fig:results-success2}
\end{figure}

\newlength{\widthfigurefailure}
\setlength{\widthfigurefailure}{0.084\linewidth}

\begin{figure}[!htbp]
\captionsetup[subfigure]{labelformat=empty}
\centering

\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/simoncelli_119.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/simoncelli_119_out_w20_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/tissus6.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/tissus6_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/simoncelli_1.png}}
\quad
\subfloat[$w_p=10, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/patch_size/simoncelli_1_out_w10_r2.png}}

\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/simoncelli_119_map_in.png}}
\quad
\subfloat[$w_p=20, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/simoncelli_119_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/tissus6_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/tissus6_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/simoncelli_1_map_in.png}}
\quad
\subfloat[$w_p=10, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/patch_size/simoncelli_1_w10_r2_map_out.png}}

\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/bark1004.png}}
\quad
\subfloat[$w_p=80, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/bark1004_out_w80_r4.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/bark1005.png}}
\quad
\subfloat[$w_p=80, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/bark1005_out_w80_r4.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/simoncelli_99.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/simoncelli_99_out_w40_r2.png}}

\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/bark1004_map_in.png}}
\quad
\subfloat[$w_p=80, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/bark1004_w80_r4_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/bark1005_map_in.png}}
\quad
\subfloat[$w_p=80, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/bark1005_w80_r4_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/simoncelli_99_map_in.png}}
\quad
\subfloat[$w_p=40, r=2$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/simoncelli_99_w40_r2_map_out.png}}

\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/IMG_9962.png}}
\quad
\subfloat[$w_p=20, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/IMG_9962_out_w20_r4.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/floor1014.png}}
\quad
\subfloat[$w_p=20, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/floor1014_out_w20_r4.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/visTex_17.png}}
\quad
\subfloat[$w_p=40, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/visTex_17_out_w40_r4.png}}

\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/IMG_9962_map_in.png}}
\quad
\subfloat[$w_p=20, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/IMG_9962_w20_r4_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/floor1014_map_in.png}}
\quad
\subfloat[$w_p=20, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/floor1014_w20_r4_map_out.png}}
\quad
\subfloat{\includegraphics[width=\widthfigurefailure]{figs/input/visTex_17_map_in.png}}
\quad
\subfloat[$w_p=40, r=4$]{\includegraphics[width=2\widthfigurefailure]{figs/results_25/output/failures/visTex_17_w40_r4_map_out.png}}

\caption{Failures of Efros and Freeman image quilting algorithms. The small images represent the example texture and the big ones the corresponding synthesis result. The examples in the first row show failures related to an incorrect patch size. The examples in the second row show the case of verbatim copy failures. The examples in the the third row show the case of growing garbage. For all the examples the patch size $w_p$ and the ratio $r$ used is indicated.}
\label{fig:results-failures}
\end{figure}

\subsubsection*{Influence of the patch size} First of all the patch size influence is analyzed and these results are shown in Figure~\ref{fig:patch-size}. For this, $w_o = 0.25$ and $\varepsilon = 0.1$ are fixed and the synthesis results for $w_p \in \lbrace 10, 20, 40, 80\rbrace $ are compared. This parameter clearly depends on the input sample. In general, macro textures have details at different scales. If the patch size chosen is not able to capture them all then the synthesis fails. A second observation is that the larger the patch size is the larger the verbatim copy zones are. For smaller patches this effect is reduced, since more ``similar'' patches are available in the patch search step and then the set of candidate patches is larger, allowing more variation.

\begin{figure}[!htbp]
\captionsetup[subfigure]{labelformat=empty,position=top}
\centering

\subfloat[input]{\includegraphics[width=0.1\linewidth]{figs/input/IMG_9898.png}}
\quad
\subfloat[$w_p=10$]{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_out_w10_r2.png}}
\quad
\subfloat[$w_p=20$]{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_out_w20_r2.png}}
\quad
\subfloat[$w_p=40$]{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_out_w40_r2.png}}
\quad
\subfloat[$w_p=80$]{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_out_w80_r2.png}}

\subfloat{\includegraphics[width=0.1\linewidth]{figs/input/IMG_9898_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_w10_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/IMG_9898_w80_r2_map_out.png}}

\subfloat{\includegraphics[width=0.1\linewidth]{figs/input/simoncelli_63.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_out_w10_r2.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_out_w20_r2.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_out_w80_r2.png}}

\subfloat{\includegraphics[width=0.1\linewidth]{figs/input/simoncelli_63_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_w10_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/simoncelli_63_w80_r2_map_out.png}}

\subfloat{\includegraphics[width=0.1\linewidth]{figs/input/tissus12.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_out_w10_r2.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_out_w20_r2.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_out_w80_r2.png}}

\subfloat{\includegraphics[width=0.1\linewidth]{figs/input/tissus12_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_w10_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_w20_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.2\linewidth]{figs/patch_size/tissus12_w80_r2_map_out.png}}

\caption{Patch size influence. Each pair 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 the examples $w_o = 0.25$ and $\varepsilon = 0.1$.}
\label{fig:patch-size}
\end{figure}

\subsubsection*{Influence of the overlap size} The overlap size is an important parameter in this method. To analyze its influence $w_p=40$ and $\varepsilon=0.1$ are fixed and the results for $w_o\in \lbrace 0.10, 0.25, 0.5 \rbrace$ are compared. The general conclusion is that larger overlap sizes often give a seamless transition between patches but are more prone to growing garbage for some texture samples as can be seen in the two first examples in Figure~\ref{fig:overlap-size} 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 at the cost of decreasing the visual quality of the results as can be seen in the third example in Figure~\ref{fig:overlap-size}.

\begin{figure}[!htbp]
\captionsetup[subfigure]{labelformat=empty,position=top}
\centering

\subfloat[input]{\includegraphics[width=0.12\linewidth]{figs/input/vegetation1006.png}}
\quad
\subfloat[$w_o=0.10$]{\includegraphics[width=0.24\linewidth]{figs/results_10/vegetation1006_out_w40_r2.png}}
\quad
\subfloat[$w_o=0.25$]{\includegraphics[width=0.24\linewidth]{figs/results_25/output/failures/vegetation1006_out_w40_r2.png}}
\quad
\subfloat[$w_o=0.50$]{\includegraphics[width=0.24\linewidth]{figs/results_50/output/vegetation1006_out_w40_r2.png}}

\subfloat{\includegraphics[width=0.12\linewidth]{figs/input/vegetation1006_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_10/vegetation1006_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_25/output/failures/vegetation1006_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_50/output/vegetation1006_w40_r2_map_out.png}}

\subfloat{\includegraphics[width=0.12\linewidth]{figs/input/vegetation1013.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_10/vegetation1013_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_25/output/failures/vegetation1013_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_50/output/vegetation1013_out_w40_r2.png}}

\subfloat{\includegraphics[width=0.12\linewidth]{figs/input/vegetation1013_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_10/vegetation1013_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_25/output/failures/vegetation1013_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_50/output/vegetation1013_w40_r2_map_out.png}}

\subfloat{\includegraphics[width=0.12\linewidth]{figs/input/visTex_4.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_10/visTex_4_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_25/output/visTex_4_out_w40_r2.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_50/output/visTex_4_out_w40_r2.png}}

\subfloat{\includegraphics[width=0.12\linewidth]{figs/input/visTex_4_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_10/visTex_4_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_25/output/visTex_4_w40_r2_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.24\linewidth]{figs/results_50/output/visTex_4_w40_r2_map_out.png}}

\caption{Overlap size influence. Each pair of rows shows from left to right the input, the synthesis results for $w_o = 0.1, 0.25, 0.5$ and the corresponding color and synthesis maps. For all the examples $w_p = 40$ and $\varepsilon = 0.1$.}
\label{fig:overlap-size}
\end{figure}

\subsubsection*{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 $\varepsilon \in \lbrace 0.05, 0.1, 0.3, 0.5, 0.7 \rbrace$ are compared. Increasing $\varepsilon$ 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~\ref{fig:tolerance} where for $\varepsilon \in \lbrace 0.5, 0.7\rbrace$
 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 $\varepsilon$ lead to very large verbatim copy areas.

\begin{figure}[!htbp]
\captionsetup[subfigure]{labelformat=empty,position=top}
\centering

\subfloat[input]{\includegraphics[width=0.08\linewidth]{figs/input/simoncelli_96.png}}
\quad
\subfloat[$\varepsilon = 0.05$]{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_out_w40_r2_t5.png}}
\quad
\subfloat[$\varepsilon = 0.1$]{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_out_w40_r2_t10.png}}
\quad
\subfloat[$\varepsilon = 0.3$]{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_out_w40_r2_t30.png}}
\quad
\subfloat[$\varepsilon = 0.5$]{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_out_w40_r2_t50.png}}
\quad
\subfloat[$\varepsilon = 0.7$]{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_out_w40_r2_t70.png}}

\subfloat{\includegraphics[width=0.08\linewidth]{figs/input/simoncelli_96_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_w40_r2_t5_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_w40_r2_t10_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_w40_r2_t30_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_w40_r2_t50_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_96_w40_r2_t70_map_out.png}}

\subfloat{\includegraphics[width=0.08\linewidth]{figs/input/simoncelli_7.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_out_w40_r2_t5.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_out_w40_r2_t10.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_out_w40_r2_t30.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_out_w40_r2_t50.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_out_w40_r2_t70.png}}

\subfloat{\includegraphics[width=0.08\linewidth]{figs/input/simoncelli_7_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_w40_r2_t5_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_w40_r2_t10_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_w40_r2_t30_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_w40_r2_t50_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/simoncelli_7_w40_r2_t70_map_out.png}}

\subfloat{\includegraphics[width=0.08\linewidth]{figs/input/wall1003.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_out_w40_r2_t5.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_out_w40_r2_t10.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_out_w40_r2_t30.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_out_w40_r2_t50.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_out_w40_r2_t70.png}}

\subfloat{\includegraphics[width=0.08\linewidth]{figs/input/wall1003_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_w40_r2_t5_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_w40_r2_t10_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_w40_r2_t30_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_w40_r2_t50_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/wall1003_w40_r2_t70_map_out.png}}

\subfloat{\includegraphics[width=0.08\linewidth]{figs/input/IMG_9873.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_out_w40_r2_t5.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_out_w40_r2_t10.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_out_w40_r2_t30.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_out_w40_r2_t50.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_out_w40_r2_t70.png}}

\subfloat{\includegraphics[width=0.08\linewidth]{figs/input/IMG_9873_map_in.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_w40_r2_t5_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_w40_r2_t10_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_w40_r2_t30_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_w40_r2_t50_map_out.png}}
\quad
\subfloat{\includegraphics[width=0.16\linewidth]{figs/tolerance/IMG_9873_w40_r2_t70_map_out.png}}

\caption{Tolerance error influence. Each couple of rows shows from left to right the input, the synthesis results for $\varepsilon = 0.05, 0.1, 0.3, 0.5, 0.7$ and the corresponding color and synthesis maps. For all the examples $w_p = 40$ and $w_o = 0.25$.}
\label{fig:tolerance}
\end{figure}

\section{Conclusion}

In this paper we analyzed in detail Efros and Freeman's texture synthesis algorithm~\cite{Efros_Freeman_image_quilting_2001}. 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 at 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~\cite{Efros_Leung_texture_synthesis_non_parametric_sampling_1999}. Another issue that arises in~\cite{Efros_Freeman_image_quilting_2001} 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 multi-core processors.
%

\section*{Image Credits}
\begin{small}

Images by the authors except:\\
\includegraphics[width=2em]{figs/input/simoncelli_8.png}
\includegraphics[width=2em]{figs/input/simoncelli_11.png}
\includegraphics[width=2em]{figs/input/simoncelli_103.png}
\includegraphics[width=2em]{figs/input/simoncelli_73.png}
\includegraphics[width=2em]{figs/input/simoncelli_29.png}
\includegraphics[width=2em]{figs/input/simoncelli_119.png}
\includegraphics[width=2em]{figs/input/simoncelli_1.png}
\includegraphics[width=2em]{figs/input/simoncelli_99.png}
\includegraphics[width=2em]{figs/input/simoncelli_63.png}
\includegraphics[width=2em]{figs/input/simoncelli_96.png}
\includegraphics[width=2em]{figs/input/simoncelli_7.png}
By Eero P. Simoncelli (\url{http://www.cns.nyu.edu/~lcv/texture/}).\\
\includegraphics[width=2em]{figs/input/IMG_9955.png}
\includegraphics[width=2em]{figs/input/IMG_9946.png}
\includegraphics[width=2em]{figs/input/IMG_9973.png}
\includegraphics[width=2em]{figs/input/IMG_9962.png}
\includegraphics[width=2em]{figs/input/IMG_9898.png}
\includegraphics[width=2em]{figs/input/IMG_9873.png}
By Yann Gousseau.\\
\includegraphics[width=2em]{figs/input/tissus4.png}
\includegraphics[width=2em]{figs/input/tissus10.png}
\includegraphics[width=2em]{figs/input/tissus6.png}
\includegraphics[width=2em]{figs/input/tissus12.png}
Unknown source, available on the internet.\\
\includegraphics[width=2em]{figs/input/visTex_2.png}
\includegraphics[width=2em]{figs/input/visTex_15.png}
\includegraphics[width=2em]{figs/input/visTex_4.png}
Fabric.0013, Food.0001 and Tile.0007 from the VisTex database (\url{http://vismod.media.mit.edu/vismod/imagery/VisionTexture/vistex.html}).\\

\end{small}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{siam}
\bibliography{biblio}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
%-------------------------------------------------------------------------------
