{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Multiplicative Cascade Synthesis of Signals and Textures", "========================================================", "", "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) 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}$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import graphics_4_multiplicative_cascade as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["1D Multiplicative Cascades", "--------------------------", "We consider here the synthesis of 1D signals using multiplicative", "compound Poisson cascades (CPC).", "", "", "Duration of the simulation : [0, T]"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["T = 10"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Resolution of the cascade (the smallest scale at which details will be added)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["rmin = 0.02"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Sampling period of the simulation. Number of points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Delta_t = rmin/ 2", "n = T/ Delta_t + 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Number of multipliers Wi."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["lambda = (1/ rmin-1)*(T + 1); % to ensure scale invariance, scales are distributed as 1/ r^2.", "N = round(lambda)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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]."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Scales ri of the points.", "Should be distributed according to dr/r^2, to have scale invariance."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["umax = 1/ rmin-1", "ui = [zeros(1, length(ti_1)) rand(1, N) * umax]; % ui = 1/ ri-1", "ri = (1 + ui).^(-1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the Poisson point process in the scale-space plane."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["figure(1)", "", "plot(ti, ri, '.')", "axis([-1/ 2, T + 1/ 2, 0, 1])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sigma2 = 0.2", "mu = -sigma2/ 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Condition of non-degeneracy:"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["if -(exp(2*(sigma2 + mu))-1) <-1", "   disp('Be careful ! This cascade will degenerate as rmin - > 0 !')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Random log-normal multipliers."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["N = length(ti);   % the number of multipliers = number of time-scale points.", "Wi = exp(randn(N, 1)*sqrt(sigma2) + mu)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Positions along time axis."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["t = linspace(0, T, n)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Initialize the signal and normalize the measure."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["H1 = 1 - exp(mu + sigma2/ 2)", "f = ones(n, 1) * exp(H1) / rmin^H1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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))|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["i = 1", "I = find(abs(t-ti(i)) <= ri(i)/ 2); % t belongs to a disk centered on ti(i)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform the multiplication with the random multiplier."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f(I) = f(I) * Wi(i)"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot(t, f)", "axis([0 T 0 1.1*max(f)])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Perform the cascade. Display intermediate steps."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo1()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the random measure."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["figure(1)", "", "plot(t, f)", "axis([0 T 0 1.1*max(f)])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Compute several realization for the same log-normal parameters."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo2()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 3__", "", "Compute realizations for different log-normal parameters |mu| and |sigma2|.", "Use the same distribution of points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo3()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["2D Multiplicative Cascades", "--------------------------", "To generate 2D cascade, one needs to throw points on a 3D scale space", "domain.", "", "", "Size of the image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 128"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Minimum scale, should be roughly |1/n|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["rmin = 1/ n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Maximum scale."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Xmax = 1", "Ymax = 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Number of points in the cascade."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Scale of the points.", "Should be distributed according to 1/height^3, to have scaling", "invariance."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Position of the points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["xi = -1/ 2 + rand(1, N) * (Xmax + 1)", "yi = -1/ 2 + rand(1, N) * (Ymax + 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the points in the scale-space plane."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["h = plot3(xi, yi, ri, '.')", "axis([-1/ 2, Xmax + 1/ 2, -1/ 2, Ymax + 1/ 2, 0, 1])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Parameters for the log-normal law."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sigma2 = 0.08", "mu = -sigma2/ 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Random log-normal multipliers."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Wi = exp(randn(N, 1)*sqrt(sigma2) + mu)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Position in the X/Y plane.", "We enlarge the square in order to be able to use periodic boundary", "conditons."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["x = linspace(0, Xmax, n)", "y = linspace(0, Ymax, n)", "[X, Y] = meshgrid(x, y)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Initialization and normalization of the image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["H1 = 1 - exp(mu + sigma2/ 2)", "f = ones(n)/ rmin^H1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["i = 1", "I = find((X-xi(i)).^2 + (Y-yi(i)).^2 <= ri(i)^2/ 4)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Multiplication of the image with the random multiplier."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f(I) = f(I) * Wi(i)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 4__", "", "Perform the full cascade, display intermediate steps."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo4()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the image. It corresponds to a 2D multi-fractal measure."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["To compute the final texture, we perform a *spectral integration*, which", "corresponds to a low pass filtering.", "", "", "Fourier frequency localizations."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["x = [0: n/ 2 -n/ 2 + 1: -1]", "[U, V] = meshgrid(x, x)", "S = sqrt(U.^2 + V.^2)", "S(1, 1) = 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Exponent of integration."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["alpha = .5"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Fourier domain integration."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["F = real(ifft2(fft2(f)./ S.^alpha))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 5__", "", "Compute the fractional integration for several values of alpha."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo5()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 6__", "", "Perform the  cascade for several log-normal parameters |mu| and |sigma2|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo6()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}