Multiplicative Cascade Synthesis of Signals and Textures

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This numerical tour explores multifractal signal and texture synthesis.

This tour is written by Pierre Chainais and Gabriel Peyr .

The processes we deal with belong to the family of Infinitely Divisible Cascades (IDC). Only the simulation of the subfamily of Compound Poisson Cascades (CPC) is really simple to implement in 1D, 2D or even ND. Indeed, the synthesis of CPC can be understood as the product of a random number of indicator function of balls (of a 1D segment or a 2D ball) with randomized radius and randomized amplitude.

If the distribution of the amplitudes and radii is well chosen, this leads to the synthesis of a function that is the density of a positive scale invariant measure. More precisely, this measure is multifractal.

To obtain the final measure signal/image, the simulated measure density is integrated or pseudo-integrated thanks to some scale invariant low-pass filtering in the Fourier domain.

The application of cascade for texture synthesis is detailed in

Infinitely divisible cascades to model the statistics of natural images, P. Chainais, IEEE Trans. on Pattern Analysis and Machine Intelligence, Vol. 29 no 12, Dec. 2007.

Visit the homepage of <http://www.isima.fr/~pchainai/PUB/software.html Pierre Chainais> for additional information and softwares.

In [2]:
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/graphics_4_multiplicative_cascade')

1D Multiplicative Cascades

We consider here the synthesis of 1D signals using multiplicative compound Poisson cascades (CPC).

Duration of the simulation : [0, T]

In [3]:
T = 10;

Resolution of the cascade (the smallest scale at which details will be added).

In [4]:
rmin = 0.02;

Sampling period of the simulation. Number of points.

In [5]:
Delta_t = rmin/2;
n = T/Delta_t+1;

Number of multipliers Wi.

In [6]:
lambda = (1/rmin-1)*(T+1); % to ensure scale invariance, scales are distributed as 1/r^2.
N = round(lambda);

N = poissrnd(N,1,1); % rigorously, N is a Poisson r.v. of expectation lambda

Time positions of the Poisson point process in the time/scale plane are uniformly distributed to ensure stationarity. Side effects are avoided by extending the cascade to [-1/2, 0] and [T T+1/2].

In [7]:
ti = -1/2+rand(1,N) * (T+1);
ti_1 = -1/2+rand(1,round(T+1))*(T+1);
ti = [ti_1 ti];   % for exact scale invariance

Scales ri of the points. Should be distributed according to dr/r^2, to have scale invariance.

In [8]:
umax = 1/rmin-1;
ui = [zeros(1,length(ti_1)) rand(1,N) * umax ]; % ui = 1/ri-1
ri = (1+ui).^(-1);

Display the Poisson point process in the scale-space plane.

In [9]:
figure(1)
clf
plot(ti, ri, '.')
axis([-1/2,T+1/2,0,1]);

Parameters for the law of multipliers Wi (Wi>0). Here we choose a log-normal law. Another simple possible choice is to set Wi=2/3 for all the Wi.

In [10]:
sigma2 = 0.2;
mu = -sigma2/2;

Condition of non-degeneracy:

In [11]:
if -(exp(2*(sigma2+mu))-1)<-1
   disp('Be careful ! This cascade will degenerate as rmin -> 0 !')
end

Random log-normal multipliers.

In [12]:
N = length(ti);   % the number of multipliers = number of time-scale points.
Wi = exp( randn(N,1)*sqrt(sigma2)+mu );

Positions along time axis.

In [13]:
t = linspace(0,T,n);

Initialize the signal and normalize the measure.

In [14]:
H1 = 1 - exp(mu+sigma2/2);
f = ones(n,1) * exp(H1) / rmin^H1;

We show here the first step of the multiplicative cascade: iterations are on the multipliers (ti,ri,Wi).

Select the points in the cone of influence of |(ti(1),ri(1))|.

In [15]:
i = 1;
I = find(abs(t-ti(i))<=ri(i)/2); % t belongs to a disk centered on ti(i)

Perform the multiplication with the random multiplier.

In [16]:
f(I) = f(I) * Wi(i);
In [17]:
clf
plot(t,f)
axis([0 T 0 1.1*max(f)])

Exercise 1

Perform the cascade. Display intermediate steps.

In [18]:
exo1()
In [19]:
%% Insert your code here.

Display the random measure.

In [20]:
figure(1)
clf
plot(t,f)
axis([0 T 0 1.1*max(f)])

Exercise 2

Compute several realization for the same log-normal parameters.

In [21]:
exo2()
In [22]:
%% Insert your code here.

Exercise 3

Compute realizations for different log-normal parameters |mu| and |sigma2|. Use the same distribution of points.

In [23]:
exo3()
In [24]:
%% Insert your code here.

2D Multiplicative Cascades

To generate 2D cascade, one needs to throw points on a 3D scale space domain.

Size of the image.

In [25]:
n = 128;

Minimum scale, should be roughly |1/n|.

In [26]:
rmin = 1/n;

Maximum scale.

In [27]:
Xmax = 1;
Ymax = 1;

Number of points in the cascade.

In [28]:
lambda = 2/pi*(1/rmin^2-1)*(Xmax+1)*(Ymax+1); % density g(r)dr=4/pi/r^3 dr
N = round(lambda); % should be a Poisson r.v. with expectation lambda.

Scale of the points. Should be distributed according to 1/height^3, to have scaling invariance.

In [29]:
umax = (1/rmin^2-1);          % u will be a uniform variable in [0 1/rmin^2-1]
ui = rand(1,N) * umax;
ri = (1/rmin^2-ui).^(-1/2);

Position of the points.

In [30]:
xi = -1/2 + rand(1,N) * (Xmax+1);
yi = -1/2 + rand(1,N) * (Ymax+1);

Display the points in the scale-space plane.

In [31]:
clf
h = plot3(xi, yi, ri, '.');
axis([-1/2,Xmax+1/2,-1/2,Ymax+1/2,0,1]);

Parameters for the log-normal law.

In [32]:
sigma2 = 0.08;
mu = -sigma2/2;

Random log-normal multipliers.

In [33]:
Wi = exp( randn(N,1)*sqrt(sigma2)+mu );

Position in the X/Y plane. We enlarge the square in order to be able to use periodic boundary conditons.

In [34]:
x = linspace(0,Xmax,n);
y = linspace(0,Ymax,n);
[X,Y]= meshgrid(x,y);

Initialization and normalization of the image.

In [35]:
H1 = 1 - exp(mu+sigma2/2);
f = ones(n)/rmin^H1;

We give here the example of the first mutiplication.

Localization of the signal locations that are influenced by the scale/space point indexed by |(xi(1),yi(1),ri(1))|. This corresponds to the intersection of the image plane and a cone of influence.

In [36]:
i = 1;
I = find( (X-xi(i)).^2+(Y-yi(i)).^2 <=ri(i)^2/4 );

Multiplication of the image with the random multiplier.

In [37]:
f(I) = f(I) * Wi(i);

Exercise 4

Perform the full cascade, display intermediate steps.

In [38]:
exo4()
In [39]:
%% Insert your code here.

Display the image. It corresponds to a 2D multi-fractal measure.

In [40]:
clf
imageplot(f)

To compute the final texture, we perform a spectral integration, which corresponds to a low pass filtering.

Fourier frequency localizations.

In [41]:
x = [0:n/2 -n/2+1:-1];
[U,V] = meshgrid(x,x); 
S = sqrt(U.^2+V.^2); 
S(1,1) = 1;

Exponent of integration.

In [42]:
alpha = .5;

Fourier domain integration.

In [43]:
F = real( ifft2( fft2(f)./S.^alpha ) );

Exercise 5

Compute the fractional integration for several values of alpha.

In [44]:
exo5()
In [45]:
%% Insert your code here.

Exercise 6

Perform the cascade for several log-normal parameters |mu| and |sigma2|.

In [46]:
exo6()
In [47]:
%% Insert your code here.