
%   This script file can be used to reproduce the experiments
%   of the paper
%     A Texton for Fast and Flexible Gaussian Texture Synthesis. 
%       (Bruno Galerne, Arthur Leclaire, Lionel Moisan), 
%        preprint MAP5 n°2014-08.
%   If you use this code for further work, please be sure to make
%   reference to this paper.


%% FIG 1 : SPOT NOISE SYNTHESIS AT LOW INTENSITY

close all
clear all

% Example of the article
u = im2double(imread('textures/RockSediment0006_7_thumbhuge_128.png'));

[M,N,~] = size(u);
mu = repmat(mean(mean(u),2),[M N 1]);
tu = 1/sqrt(M*N)*(u-mu);

% Creation of the support
r = 15; % support half-size
M2 = floor(M/2)+1; N2 = floor(N/2)+1;
sup = zeros(M,N);
xind = M2-r:M2+r;
yind = N2-r:N2+r;
sup(xind,yind) = Inf;
suprgb = repmat(sup,[1 1 3]);
cardsup = sum(sum(sup>0));

% SO texton computation
options.err = 0;
options.nmax = 1000;
options.visu = 1;
sot = so_texton(u,suprgb,options);
% Crop the central part
sotc = sot(xind,yind,:);
[m n C] = size(sotc);

% Synthesis size
Ms = 128; Ns = 128;
% Mean image for synthesis
mus = repmat(mu(1,1,:),[Ms Ns 1]);

% DSN
[sn30,pp30] = dsn(sotc,mus,30/cardsup);
% % ADSNs
% adsnu = adsn(tu,mus);
% adsnsot = adsn(sotc,mus);

% Normalization of the SOT for visualization
sotvisu = mu + sqrt(cardsup)*sot;
% cropped version : 
% sotvisu = mu(1:m,1:n,:)+sqrt(cardsup)*sotc;

% Visualization
figure
subplot(1,4,1)
imshow(u)
title('Original texture u')
subplot(1,4,4)
imshow(sn30)
title('SN(t), 30 impacts/pix')
subplot(1,4,3)
imshow(pp30,[])
title('Underlying Poisson Process')
subplot(1,4,2)
imshow(sotvisu)
title(['SOT t of size ' num2str([m n])])

%% FIG 2 : ITERATES OF ALGORITHM 2

close all
clear all

% Decomment and execute the following lines to run 
%       the Monte-Carlo simulation

% % Initialize random seed
% RandStream.setDefaultStream(RandStream('mt19937ar','seed',1));
% % Original image
% u = im2double(rgb2gray(imread('textures/IMG_7396.png')));
% u = u(1:256,1:384);
% [M,N] = size(u);
% figure
% imshow(u)
% title('Original texture')
% % Parameters
% niter = 1000;
% options.err = 0;
% options.nmax = niter;
% % Choice of a circular support
% renv = 15;
% M2 = floor(M/2)+1; N2 = floor(N/2)+1;
% sup = disc(M,N,M2,N2,renv,1);
% sup(sup>0) = Inf;
% % Computation of several SOTs with RPN initializations
% nb = 100;
% dot = zeros(nb,niter);
% err = zeros(nb,niter);
% for k=1:nb
%     disp(num2str(k))
%     [~,evodot,evoerr] = so_texton(u,sup,options);
%     dot(k,:) = evodot;
%     err(k,:) = evoerr;
% end

% Or load the following file containing the data
load vary_tex6_ostexton_iter1000.mat

% Save input image :
mdot = mean(dot);
sdot = sqrt(var(dot));
x = 1:100; 

figure;
plot(x,mdot(x),'LineWidth',2)
xlabel('n')
ylabel('Mean')

figure;
plot(x,sdot(x),'LineWidth',2)
xlabel('n')
ylabel('Standard Deviation')

%% FIG 3 : SYNTHESIS OF COLOR TEXTURES

close all
clear all

tex = 1; % texture choice
switch tex
    % Examples used in the article :
    case 1
        % Texture 1 :
        u = im2double(imread('textures/IMG_7396.png')); u = u(1:256,1:384,:);
        [M,N,~] = size(u);
        mu = repmat(mean(mean(u),2),[M N 1]);
        u = mu+2*(u-mu); % Increase the contrast
    case 2
        u = im2double(imread('textures/wall1011.png')); u = u(21:276,21:404,:);
    case 3
        u = im2double(imread('textures/c28.png')); u = u(1:256,1:384,:);
   
    % Other examples :
    case 4
        u = im2double(imread('textures/c05.png')); u = u(1:256,1:384,:);
    case 5
        u = im2double(imread('textures/c17.png')); u = u(1:256,1:384,:);
    case 6
        u = im2double(imread('textures/c10.png')); u = u(1:256,1:384,:);
end

[M,N,~] = size(u);
mu = repmat(mean(mean(u),2),[M N 1]);
tu = 1/sqrt(M*N)*(u-mu);

% Replace u by its RPN. 
u = rpn(tu,mu);

figure
imshow(u)
title('Original random phase texture')

% Creation of the support
r = 15; % support half-size (15 for tex 1 and 25 for the others).
M2 = floor(M/2)+1; N2 = floor(N/2)+1;
sup = zeros(M,N);
xind = M2-r:M2+r;
yind = N2-r:N2+r;
sup(xind,yind) = Inf;
suprgb = repmat(sup,[1 1 3]);
cardsup = sum(sum(sup>0));

% SO texton computation
options.err = 1e-3;
options.nmax = 1000;
options.visu = 1;
sot = so_texton(u,suprgb,options);
sotc = sot(xind,yind,:); % Crop the central part
[m n C] = size(sotc);

% Covariance correction
sotc_cc = cov_correction(sotc,u);

% Visual comparison of the textons
muc = mu(xind,yind,:);
figure
subplot(1,2,1)
imshow(muc+sqrt(m*n)*sotc)
title(sprintf(['SOT t \n RME = ' ...
    num2str(otdist(tu,sotc,M,N,1))]))
subplot(1,2,2)
imshow(muc+sqrt(m*n)*sotc_cc)
title(sprintf(['SOT t_{cc} with color correction \n RME = ' ...
    num2str(otdist(tu,sotc_cc,M,N,1))]))

% Synthesis
Ms = M; 
Ns = N;
mus = repmat(mu(1,1,:),[Ms Ns 1]);

% DSN
nimp = 50; % number of impacts per pixel
sn = dsn(sotc,mus,nimp/cardsup);
sn_cc = dsn(sotc_cc,mus,nimp/cardsup);

% ADSNs
adsnsot = adsn(sotc,mus);
adsnsot_cc = adsn(sotc_cc,mus);

% Is the covariance of a pixel marginal ok with color correction ?
cov(reshape(u,[M*N 3]))
cov(reshape(adsnsot_cc,[M*N 3]))
% yes
    
% Visu
figure
imshow(u)
title('Random phase texture u')
figure
imshow(adsnsot)
title('ADSN(t)')
figure
imshow(sn)
title(['SN(t), ' num2str(nimp) ' impacts/pix'])
figure
imshow(adsnsot_cc)
title('ADSN(t_{cc})')
figure
imshow(sn_cc)
title(['SN(t_{cc}), ' num2str(nimp) ' impacts/pix'])


%% FIG 4 : DSN SYNTHESIS OF A NATURAL COLOR TEXTURE, COMPARISON

close all
clear all

u = im2double(imread('textures/moquette_10_128x128.png'));

[M,N,~] = size(u);
mu = repmat(mean(mean(u),2),[M N 1]);
tu = 1/sqrt(M*N)*(u-mu);

% Creation of the support
r = 15; % support half-size
M2 = floor(M/2)+1; N2 = floor(N/2)+1;
sup = zeros(M,N);
xind = M2-r:M2+r;
yind = N2-r:N2+r;
sup(xind,yind) = Inf;
cardsup = sum(sum(sup>0));

% SO texton computation
suprgb = repmat(sup,[1 1 3]);
options.err = 0;
options.nmax = 1000;
options.visu = 1;
sot = so_texton(u,suprgb,options);

% Crop the central part
sotc = sot(xind,yind,:);
muc = mu(xind,yind,:);
[m n ~] = size(sotc);

% For baseline comparison
% Cropped rpn
rpnu = rpn(tu,zeros(size(u)));
trpn = 1/sqrt(m*n)*rpnu(1:m,1:n,:);
% Luminance texton
tlum = circshift(luminance_texton(u),[floor(M/2) floor(N/2)]);
tlumc = tlum(xind,yind,:);

% Visual comparison of the textons
figure
subplot(1,3,1)
imshow(muc+sqrt(m*n)*sotc)
title(sprintf(['Synthesis-Oriented Texton t \n RME = ' ...
    num2str(otdist(tu,sotc,M,N,1))]))
subplot(1,3,2)
imshow(muc+sqrt(m*n)*tlumc)
title(sprintf(['Cropped Luminance Texton t_{lum}^c \n RME = ' ...
    num2str(otdist(tu,tlumc,M,N,1))]))
subplot(1,3,3)
imshow(muc+sqrt(m*n)*trpn)
title(sprintf(['Cropped RPN Texton t_{rpn} \n RME = ' ...
    num2str(otdist(tu,trpn,M,N,1))]))
maximize

% Synthesis

Ms = 128; Ns = 128;
mus = repmat(mu(1,1,:),[Ms Ns 1]);

% SN
sn10 = dsn(sotc,mus,10/cardsup);
sn30 = dsn(sotc,mus,30/cardsup);
sntrpn = dsn(trpn,mus,30/cardsup);
sntlum = dsn(tlumc,mus,30/cardsup);
% ADSNs
adsnu = adsn(tu,mus);
adsnsot = adsn(sotc,mus);
adsntlum = adsn(tlumc,mus);

% Visualization
figure
subplot(2,4,1)
imshow(u)
title('Original texture u')
subplot(2,4,2)
imshow(sn10)
title('DSN(t), 10 impacts/pix')
subplot(2,4,3)
imshow(sn30)
title('DSN(t), 30 impacts/pix')
subplot(2,4,4)
imshow(adsnsot)
title('ADSN(t)')
subplot(2,4,5)
imshow(adsnu)
title('ADSN(t_u)')
subplot(2,4,6)
imshow(sntlum)
title('DSN(t_{lum}^c), 30 impacts/pix')
subplot(2,4,7)
imshow(sntrpn)
title('DSN(t_{rpn}), 30 impacts/pix')
maximize


%% FIG 5 : INFLUENCE OF THE SUPORT SIZE

close all
clear all

u = im2double(imread('textures/ConcreteFloors0068_2_thumbhuge_128.png'));
[M,N,C] = size(u);
M2 = floor(M/2)+1; N2 = floor(N/2)+1;
mu = repmat(mean(mean(u),2),[M N 1]);
u = mu + 2*(u-mu); % increase the contrast
un = 1/sqrt(M*N)*(u-mu);

figure
imshow(u)
title('Original texture')

% SO texton parameters
r = [5 15 25];
options.err = 0;
options.nmax = 1000;
% DSN synthesis parameter
nb = 50;

adsnu = adsn(un,mu);
tlum = circshift(luminance_texton(u),[floor(M/2) floor(N/2)]);

adsnsot_tab = cell(size(r));
adsntlumc_tab = cell(size(r));

dsot = zeros(size(r));
dlumt = zeros(size(r));

for k=1:length(r)
sup = disc(M,N,M2,N2,r(k),Inf);
suprgb = repmat(sup,[1 1 3]);
xind = max(sup,[],2)>0; m = sum(xind(:)>0);
yind = max(sup)>0; n = sum(yind(:)>0);
cardsup = sum(sum(sup>0));
% SO texton
sot = so_texton(u,suprgb,options);
sotc = sot(xind,yind,:);
% Cropped luminance texton
tlumc = tlum(xind,yind,:);
% Syntheses
snsot = dsn(sotc,mu,nb/cardsup);
adsnsot = adsn(sot,mu);
adsnsot_tab{k} = adsnsot;
sntlumc = dsn(tlumc,mu,nb/cardsup);
adsntlumc = adsn(tlumc,mu);
adsntlumc_tab{k} = adsntlumc;
    % Visualization
    figure
    subplot(2,3,1)
    imshow(u)
    title('Original')
    subplot(2,3,2)
    imshow(adsnsot)
    title(['ADSN(t^{' num2str(r(k)) '})'])
    subplot(2,3,3)
    imshow(adsntlumc)
    title(['ADSN(t_{lum}^{' num2str(r(k)) '})'])
    subplot(2,3,5)
    imshow(snsot)
    title(['SN(t^{' num2str(r(k)) '})'])
    subplot(2,3,6)
    imshow(sntlumc)
    title(['SN(t_{lum}^{' num2str(r(k)) '})'])
    drawnow
% Relative model error
dsot(k) = otdist(un,sot,M,N,1);
dlumt(k) = otdist(un,tlumc,M,N,1);
end

figure
plot(r,dsot,'+r',r,dlumt,'ob')
xlabel('r')
legend({['RME(t_u,t^{' num2str(r(k)) '})'],...
    ['RME(t_u,t_{lum}^{' num2str(r(k)) '})']});

figure
subplot(1,length(r)+1,1)
imshow(adsnu)
title('ADSN(t_u)')
for k=1:length(r)
    subplot(1,length(r)+1,1+k)
    imshow(adsnsot_tab{k})
    title(sprintf(['ADSN(t^{' num2str(r(k)) '}) \n RME = ' num2str(dsot(k))]));
end
maximize

figure
subplot(1,length(r)+1,1)
imshow(adsnu)
title('ADSN(t_u)')
for k=1:length(r)
    subplot(1,length(r)+1,1+k)
    imshow(adsntlumc_tab{k})
    title(sprintf(['ADSN(t_{lum}^{' num2str(r(k)) '}) \n RME = ' num2str(dlumt(k))]));
end
maximize
