% Correction TP séance 2
% cours CES8AH - Ecole des Mines de Nancy
% Frederic Sur 2016


% 1. Manipulations élémentaires

clear all, close all

Ima=imread('boat.tif');
% on essaiera aussi avec d'autres images, comme D3:
% Ima=imread('D3.tif');

figure, imshow(Ima);
% imshow montre les images en int

% on préfèrera convertir l'image en double
Ima=double(Ima);
figure;
colormap(gray);
imagesc(Ima); 
axis equal

fIma=fft2(Ima);
% spectre de puissance:
figure, imagesc(abs(fftshift(fIma)).^2), colormap(gray), axis equal
% on ne voit pas grand-chose

% en échelle log:
figure, imagesc(log(abs(fftshift(fIma)))), colormap(gray), axis equal
figure, imagesc(log(abs(fIma))), colormap(gray), axis equal

% vérification de la formule de Parseval:
sum(abs(Ima(:)).^2)-mean(abs(fIma(:)).^2)

% influence des discontinuités aux bords sur le spectre:

H=hann(size(Ima,1))*hann(size(Ima,2))';
figure, imagesc(H), colormap(gray), axis equal, colorbar
figure, imagesc(H.*Ima), colormap(gray), axis equal
fIma2=fft2(Ima.*H);
figure, imagesc(log(abs(fftshift(fIma2)))), colormap(gray),

% Hautes fréquences et basses fréquences
halfwidth=50;  % (essayer avec différentes valeurs)
G=zeros(size(Ima));
for i=1+size(Ima,1)/2+[-halfwidth:halfwidth],
  for j=1+size(Ima,2)/2+[-halfwidth:halfwidth],
    G(i,j)=1;
  end; 
end;

figure, imagesc(log(abs(G.*fftshift(fIma)))), colormap(gray), axis equal
figure, imagesc(log(abs((1-G).*fftshift(fIma)))), colormap(gray), axis equal

ImaBF=ifft2(ifftshift(G.*fftshift(fIma)));
ImaHF=ifft2(ifftshift((1-G).*fftshift(fIma)));

figure, imagesc(ImaBF), colormap(gray), colorbar
figure, imagesc(ImaHF), colormap(gray), colorbar




% 2. Etudes d'images de textures

clear all, close all

% Expérience d'Openheim et Lim

Ima1=imread('boat.tif');
Ima2=imread('house.tif');

figure; imshow(Ima1); title('Ima1');

figure; imshow(Ima2); title('Ima2');

fIma1=fft2(Ima1);
fIma2=fft2(Ima2);

Ima_m1_p2 = ifft2( abs(fIma1).*exp(i*angle(fIma2)) ,'symmetric');
Ima_m2_p1 = ifft2( abs(fIma2).*exp(i*angle(fIma1)) ,'symmetric');

figure; imshow(uint8(Ima_m1_p2)); title('module Ima1, phase Ima2');
figure; imshow(uint8(Ima_m2_p1)); title('module Ima2, phase Ima1');


% synthèse de texture

v=1; % (1: mode "verbose" dans randomize_phase)

% choisir une des images (exécutez l'UNE des lignes suivantes), 
% puis lancer la suite du code
close all, Ima=double(imread('D3.tif'));
close all, Ima=double(imread('D21.tif'));
close all, Ima=double(imread('D82.tif'));
close all, Ima=double(imread('D95.tif'));
close all, Ima=double(imread('D104.tif'));
close all, Ima=double(imread('woodpost.png')); Ima=mean(Ima,3);
close all, Ima=double(imread('concretewall.png')); Ima=mean(Ima,3);
close all, Ima=double(imread('wood.png')); Ima=mean(Ima,3);
close all, Ima=double(imread('paperwall.png')); Ima=mean(Ima,3);

figure, imshow(uint8(Ima)), title('image originale');
Ima_alea=randomize_phase(Ima,v);
figure, imshow(Ima_alea); title('texture synthétisée');

% en couleur (choisir une des images)
close all, Ima=double(imread('woodpost.png'));
close all, Ima=double(imread('concretewall.png')); 
close all, Ima=double(imread('wood.png')); 
close all, Ima=double(imread('paperwall.png'));

figure, imshow(uint8(Ima)), title('texture originale');
Ima_alea=randomize_phase(Ima,v);
figure, imshow(Ima_alea); title('texture synthétisée');


% méthode d'extension de l'image:
% cf: B. Galerne, Y. Gousseau and J.-M. Morel, 
% Random Phase Textures: Theory and Synthesis, 
% IEEE Trans. on Image Processing, 2011.  

% création d'une "smooth transition function"
alpha=0.05;
WX=ones(size(Ima,1),1);
WX(1:floor(alpha*size(Ima,1)))=0;
WX(end-floor(alpha*size(Ima,1)):end)=0;
h=gausswin(floor(3*alpha*size(Ima,1)));
WX=conv(WX,h,'same');
WY=ones(size(Ima,2),1);
WY(1:floor(alpha*size(Ima,2)))=0;
WY(end-floor(alpha*size(Ima,2)):end)=0;
h=gausswin(floor(3*alpha*size(Ima,2)));
WY=conv(WY,h,'same');
W=WX*WY';
W=W/max(W(:));
W=repmat(W,[1,1,size(Ima,3)]);

m=mean(mean(Ima,1),2);
newsize=[1000,1000]; % à modifier
Imab=sqrt(prod(newsize)/(size(Ima,1)*size(Ima,2)))*W.*(Ima-repmat(m,size(Ima,1),size(Ima,2),1));
Ima2=repmat(m,newsize(1),newsize(2),1)+padarray(Imab,floor((newsize-[size(Ima,1),size(Ima,2)])/2));

figure, imshow(uint8(Ima2)), title('texture originale')
Ima_alea=randomize_phase(Ima2,v);
figure, imshow(Ima_alea); title('texture synthétisée')

