Synthesizing new images using Gaussian models is probably the most basic idea in image synthesis. However, even if it is simple, this framework is already very rich and allows many nice recipes, going from texture synthesis to texture inpainting and texture mixing. The general model is simple and explicit, but has also many subtleties, when dealing for instance with the image boundaries or color images.
This practical session shows how to obtain texton and Gaussian synthesis of color and gray textures, as well as texture mixing.
Authors:
Below is a list of packages needed to implement the algorithms.
Numpy
PIL.Image, matplotlib.pyplot
(display of images and graphics)os
(interactions with the operating system)widgets
(display result of Gaussian texture mixing)from PIL import Image
import matplotlib.pyplot as plt
import os
import numpy as np
from google.colab import widgets
To import the solutions, execute the following cell. If you are using a Windows system, comment the os.system
line, download the file by hand, and place it in the same folder as the notebook.
#os.system("wget -nc https://raw.githubusercontent.com/storimaging/Notebooks/main/ImageGeneration/Solutions/GausianTSandMixing.py")
#from GausianTSandMixing import *
Import the helper functions for loading and displaying images by running the cell below.
os.system("wget -nc https://raw.githubusercontent.com/storimaging/Notebooks/main/ImageGeneration/AuxiliarFunctions/gaussianTSandMixing.py")
from gaussianTSandMixing import *
In the next section we will upload texture images. Not all texture images are well represented by Gaussian models. Only the so-called micro-textures are. When the geometric content is strongly constrained, the Gaussian model is not able to capture it.
Try running the exercises with micro-textures, but also with regular textures like grapes.jpg or tannat.jpg to see the difference.
texture_imgnames = ["bois.png", "mur.png", "tissu.png", "nuages.png","wall1003.png", "tulum.jpg", "grapes.jpg", "tannat.jpg"]
for fname in texture_imgnames:
os.system("wget -c https://raw.githubusercontent.com/storimaging/Images/main/Textures/"+fname)
img = Image.open(fname)
In this cell we choose the texture for the next exercise. We can choose from the images uploaded in the previous cell, by changing the parameter of the loadImage
function. Once we have the image, we will also generate its grayscale version.
u = loadImage("nuages.png")
# Convert image to gray-scale
u_gray = 0.33*u[:,:,0] + 0.5*u[:,:,1] + 0.17*u[:,:,2]
# Plot both images
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))
axes[0].imshow(u)
axes[0].set_title('image')
axes[0].axis('off')
axes[1].imshow(u_gray,cmap='gray',vmin=0,vmax=1)
axes[1].set_title('gray version of image')
axes[1].axis('off')
fig.tight_layout()
plt.figure()
<Figure size 640x480 with 0 Axes>
<Figure size 640x480 with 0 Axes>
The Fourier transform of an image is very sensitive to its discontinuities. By applying the transformation, it makes the implicit assumption that the image is somehow “periodic”. However, an image is generally not periodic and the left/right borders (resp. top/bottom) of the image don’t "coincide". As a result, a "cross" appears in the image of the amplitude of the Fourier transform.
To remove this artifact, Lionel Moisan proposed in [1] to consider the decomposition of an image $u$ into $u = p + s$, where $p$ is the periodic component and $s$ is the smooth component. $p$ will have almost the same Fourier transform as $u$, without the cross effect due to the non-periodic boundaries of $u$.
The periodic component $p$ of $u$ defined over a domain $\Omega$ is defined to be the unique solution of:
$$ \begin{cases} \Delta p = \Delta_i u \\ mean(p) = mean(u) \end{cases} $$The Laplacian $\Delta$ of an image $f$ defined over a domain $\Omega$ can be defined as:
$$ \Delta f(x) = 4f(x) - \sum_{y \in N_x} f(y) $$where $N_x$ denote the neighborhood of $x \in \Omega$ for 4-connectivity.
The Laplacian $\Delta_i$ (here index $i$ means “interior”) of an image $f$ defined over a domain $\Omega$ can be defined as:
$$ \Delta_i f(x) = |N_x \cap \Omega|f(x) - \sum_{y \in N_x \cap \Omega} f(y)$$[1] Lionel Moisan. "Periodic plus smooth image decomposition." Journal of Mathematical Imaging and Vision 39.2 (2011): 161-179.
$s = u − p$ is “smooth”: its Laplacian is equal to $0$ except at the boundary. It's known that when $s$ is extended by periodicity, we have in the Fourier domain:
$$ \forall \xi \in \Omega, \; \widehat{\Delta s}(\xi) = \left( 2 \cos \left(2\pi\langle e_1,\xi \rangle \right) + 2 \cos \left( 2\pi \langle e_2, \xi \rangle \right) - 4 \right) \widehat{s}(\xi)$$Also, as $s = u − p$, we have:
$$ \Delta s = \Delta u - \Delta p = \Delta u - \Delta_i u:= B(u) $$We can compute the “boundary discontinuity” $B(u)$ using the definitions of the Laplacian and the interior Laplacian by:
$$ \forall x \in \Omega,\; B(u)(x) = \sum_{y \in N_x \cap \Omega^c} \left( u(x) -u(y) \right)$$where $\Omega^c$ denote the complementary set of $\Omega \in \mathbb{Z}^2$ and $u$ is extended to $\mathbb{Z}^2$ by periodicity.
Therefore:
$$ \widehat{\Delta s} = \widehat{B}(u) = \left( 2 \cos \left(2\pi\langle e_1,\xi \rangle \right) + 2 \cos \left( 2\pi \langle e_2,\xi \rangle \right) - 4 \right) \widehat{s}(\xi) $$and then:
$$\widehat{s}(\xi) = \frac{\widehat{B(u)}(\xi)}{ 2 \cos \left(2\pi\langle e_1,\xi \rangle \right) + 2 \cos \left( 2\pi \langle e_2,\xi \rangle \right) - 4 } $$Only the mean value of the image can not be recovered from its Laplacian, so $\widehat{s}(0) = 0$. Finally the periodic component can be computed as $p = u − s$.
Excercise 1. Write a function ppluss
that takes as input a gray-scale image $u$ and then computes its periodic and smooth components.
P,S = ppluss(u_gray)
Below we show the original image, the perioric component and the smooth component, and their corresponding Fourier transform.
printImages(u_gray, P, S, 20, 20)
printFTImages(u_gray, P, S, 20, 20)
<Figure size 640x480 with 0 Axes>
<Figure size 640x480 with 0 Axes>
In this section we will compute the texton of a gray-scale texture and then generate a Gaussian texture synthesis from it. For the exercise we will use the periodic component of the selected original texture, therefore we must start with the first task, to fulfill all the prerequisites.
Excercise 1. Write a function ChangeToPeriodic
that takes as input gray-scale image $u_0$ of size $n_h × n_w$. Replace $u_0$ with its periodic component. Compute its empirical mean $m_{u_0}$ and remove $m_{u_0}$ from the image to have a zero-mean image $u_0$. Return the new image and $m_{u_0}$ (we must keep $m_{u_0}$ to add it to the Gaussian Synthesis at the end).
You can verify that the mean of the original image and of its periodic component are the same, so the order to compute the mean doesn't affect the results.
# Load image and convert it to gray-scale
u = loadImage("mur.png")
u_gray = 0.33*u[:,:,0] + 0.5*u[:,:,1] + 0.17*u[:,:,2]
u_gray_processed, mean_u_gray = ChangeToPeriodic(u_gray)
Next we plot the histograms of the original gray-scale image and of the processed image obtained after application of the ChangeToPeriodic
function (periodic component of the original image with mean removed).
# Histogram of u_gray
Y = np.reshape(u_gray, u_gray.shape[0]*u_gray.shape[1])
plt.hist(Y,256)
plt.title('Histogram of u_gray')
plt.show()
# Histogram of u_gray_processed
X = np.reshape(u_gray_processed, u_gray_processed.shape[0]*u_gray_processed.shape[1])
plt.hist(X,256)
plt.title('Histogram of u_gray_processed')
plt.show()
Let $u_0$ be an image defined on $\Omega$ with zero mean and periodic. The texton $T(u_0)$ of $u_0$ is the image defined on $\Omega$ such that the fourier transform of $T(u_0)$ is equal to the modulus of the fourier transform of $u_0$.
$$ \widehat{T(u_0)} = | \widehat{u_0}| $$Excercise 2. Compute and visualize the texton $T(u_0)$ of $u_0$, placing 0 in the center of the image. To better see the concentration of the texton, draw the signals $x \rightarrow T(u_0)(x,0)$ and $y \rightarrow T(u_0)(0,y)$, which are two slices of the texton image passing through the origin.
# Compute Texton
STI = ComputeTexton(u_gray_processed)
# Plot Texton placing 0 at the center of the image.
real_STI = np.real(STI)
printGrayTexton(real_STI, mean_u_gray)
# Plot Slices
PlotSlices(real_STI)
<Figure size 640x480 with 0 Axes>
Consider $u_0$ an image with zero mean defined on $\Omega$ and $W$ a white noise image of variance $σ^2 = \frac{1}{|\Omega|}$ on $\Omega$. The random image $U = u_0 ∗ W$ follows a stationary Gaussian distribution that has zero mean and the same empirical covariance as $u_0$ .
We will use the texton to compute the image synthesis. It is possible to compute a synthesis larger than the original image. To compute a synthesis $v$ in $\Omega'$ ($\Omega'$ can be the same as $\Omega$ or even larger) the steps are:
Note that the normalization step can be eliminated if the noise is generated with variance $σ^2 = \frac{1}{|\Omega|}$.
Excercise 3. Generate two Gaussian textures using the texton computed on excercise 2: the first one of size $n_h × n_w$ and the second one of size $2n_h × 2n_w$.
Gaussian_texture = GenerateGaussianTexture(STI, mean_u_gray)
Gaussian_texture_2 = GenerateGaussianTextureDouble(STI, mean_u_gray)
# Plot both Textures
PlotTextures(Gaussian_texture, Gaussian_texture_2)
In this section we will compute the texton of a color texture and then generate a Gaussian texture synthesis from it. For the exercise we will use the periodic component of the selected original texture, therefore we must start with the first task, to fulfill all the prerequisites.
Excercise 1. Take a color image $u_0$ of size $n_h × n_w$. Replace it by its periodic component (i.e. replace each channel by its periodic component). Compute its empirical mean $m_{u_0} \in \mathbb{R}^3$ and remove it from the image to have a zero-mean color image. Return the new image and $m_{u_0}$ (we must keep $m_{u_0}$ to add it to the Gaussian Synthesis at the end).
u = loadImage("mur.png")
u_processed, mean_u= ChangeToPeriodic_Color(u)
Next we plot the histograms of the first channel of the color image and the first channel of the processed image by the function ChangeToPeriodic_Color
.
# Histogram of u
Y = np.reshape(u[:,:,0], u.shape[0]*u.shape[1])
plt.hist(Y, bins= 256, range=(0,1))
plt.title('Histogram of red channel of u')
plt.show()
# Histogram of u_processed
X = np.reshape(u_processed[:,:,0], u_processed.shape[0]*u_processed.shape[1])
plt.hist(X, bins = 256)
plt.title('Histogram of red channel of u_processed')
plt.show()
When $u$ is a color image, the texton $T(u)$ can be obtained by substracting the same phase from all three channels, this phase being that of a linear combination of the channels.
More precisely, if $u = (u^R,u^G,u^B)$ is a color image defined on $\Omega$, for any $α = (\alpha^R,\alpha^G,\alpha^B) \in \mathbb{R}^3$, one can define the $\alpha$-texton of $u$ by the color image $T_\alpha(u) = (T_\alpha(u)^R,T_\alpha(u)^G,T_\alpha(u)^B) $ where for any color channel $c \in \{R,G,B\}$,
$$ \forall \xi \in \Omega, \; \widehat{T_\alpha(u)^c}(\xi) = \widehat{u^c}(\xi)e^{-i\phi_{\alpha.u}(\xi)} $$where $\phi_{\alpha.u}(\xi)$ is the phase of $\widehat{\alpha.u}(\xi)$, and the image $\alpha.u = \alpha^Ru^R + \alpha^Gu^G + \alpha^Bu^B $ being a real-valued image.
So, a way to obtain the texton is:
Excercise 2. Compute and visualize the color texton image $T_\alpha(u)$ for instance with $\alpha = (0.33, 0.5, 0.17)$ (the luminance channel), placing 0 in the center of the image.
# Compute Texton
STI_color = ComputeTexton_Color(u_processed)
# Plot Texton placing 0 at the center of the image.
real_STI_color = np.real(STI_color)
printColorTexton(real_STI_color, mean_u)
WARNING:matplotlib.image:Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Excercise 3. Generate a Gaussian color texture synthesis of size $n_h × n_w$ using the texton computed on exercise 2. To implement it, apply the same idea of Gray-scale synthesis to each channel, but use the same white noise for all.
Gaussian_texture_color = GenerateGaussianTexture_Color(STI_color, mean_u)
plt.imshow(np.clip(np.real(Gaussian_texture_color),0,1))
plt.axis('off')
(-0.5, 767.5, 511.5, -0.5)
Let $u_0$ and $u_1$ be two gray-scale images defined on the same domain $\Omega$. Let $GT(u_0)$ and $GT(u_1)$ be their respective Gaussian texture model. The 2-Wasserstein distance between $GT(u_0)$ and $GT(u_1)$ can be defined as:
$$ W^2_2(GT(u_0),GT(u_1)) = \| m_{u_0} - m_{u_1}\|^2 + \frac{1}{\Omega} \sum_{\xi \neq 0} (|\widehat{u_0}(\xi)| - |\widehat{u_1}(\xi)|)^2 \\ = |\Omega| (( m_{u_0} - m_{u_1} )^2 + \|T(u_0 - m_{u_0}) - T(u_1 - m_{u_1}) \|^2)$$where $T(u_i −m_{u_i})$ is the texton of the zero mean version of $u_i$, for i = 0 or 1.
Moreover if $\rho \in [0, 1]$, then the barycenter between $GT(u_0)$ and $GT(u_1)$ with weights $(1 − \rho, \rho)$ is given by the Gaussian texture model $GT(u_\rho)$ with
$$ u_\rho =(1−\rho)(m_{u_0} +T(u_0 −m_{u_0}))+ \rho( m_{u_1} +T(u_1 −m_{u_1})) $$Excercise 1. Take two gray-scale images $u_0$ and $u_1$ of same size $n_h × n_w$. Compute their mean, remove it and replace each image by its periodic component. Then generate a sequence of $n_\rho + 1$ (with $n_\rho$ = 5 for instance) texture images from the Wasserstein barycenter Gaussian texture model $GT(u_\rho)$, for $\rho = 0, \frac{1}{n_p} , \frac{2}{n_p} , . . . , 1$.
# Read images
u_0_gray = np.copy(u_gray)
u_1 = loadImage("tulum.jpg")
u_1_gray = 0.21*u_1[:,:,0] + 0.72*u_1[:,:,1] + 0.07*u_1[:,:,2]
# Rho
n_p = 5
# Get a list of texture mixing with differents values of rho
GTrho_list = GaussianTextureMixing_gray(u_0_gray, u_1_gray, n_p)
# Compare images
compare_images(GTrho_list, n_p, gray=True)
Let $u_0$ and $u_1$ be two color images defined on the same domain $\Omega$. Let $GT(u_0)$ and $GT(u_1)$ be their respective Gaussian texture model. The 2-Wasserstein distance between $GT(u_0)$ and $GT(u_1)$ can be defined as:
$$ W^2_2(GT(u_0),GT(u_1)) = \| m_{u_0} - m_{u_1}\|^2 + \frac{1}{\Omega} \sum_{\xi \neq 0} ( \|\widehat{u_0}(\xi)\|^2 + \|\widehat{u_1}(\xi)\|^2 -2|\widehat{u_0}(\xi)^*\widehat{u_1}(\xi)| )$$Moreover if $\rho \in [0, 1]$, then the barycenter between $GT(u_0)$ and $GT(u_1)$ with weights $(1 − \rho, \rho)$ is given by the Gaussian texture model $GT(u_\rho)$ with:
$$ m_{u_\rho} = (1 - \rho)m_{u_0} + \rho m_{u_1} \;\;\; \text{ and } \;\;\; \widehat{u_\rho}(\xi) =(1−\rho)\widehat{u_0}(\xi) + \rho \frac{\widehat{u_1}(\xi)^*\widehat{u_0}(\xi)}{|\widehat{u_1}(\xi)^*\widehat{u_0}(\xi)|}\widehat{u_1}(\xi) $$Excercise 2. From two color images, generate a sequence of $n_\rho + 1$ (with $n_\rho$ = 5 for instance) texture images from the Wasserstein barycenter Gaussian texture model $GT(u_\rho)$, for $\rho = 0, \frac{1}{n_p} , \frac{2}{n_p} , . . . , 1$.
# Read images
u_0 = np.copy(u)
u_1 = loadImage("tulum.jpg")
# Rho
n_p = 5
# Get a list of texture mixing with differents values of rho
Grho_list = GaussianTextureMixing_color(u_0, u_1, n_p)
# Compare images
compare_images(Grho_list, n_p, gray=False)