% Partie I

% EOS 350D 100ISO

clear all, close all
Ima=double(imread('mire_EOS350D_R100.png'));
figure, imagesc(Ima), colormap(gray)
moyenne=[];
variance=[];

Row={100:150,200:280,330:380};
Col={70:120,180:240,290:350,410:460,520:580};

for x=1:numel(Row),
    for y=1:numel(Col),
        bloc=Ima(Row{x},Col{y});
        moyenne=[moyenne; mean(bloc(:))];
        variance=[variance;var(bloc(:))];
    end
end

figure, plot(moyenne, variance,'x'), hold on
[coef,coefint]=regress(variance,[ones(length(variance),1),moyenne]);
X=[min(moyenne):10:max(moyenne)];
plot(X,coef(1)+coef(2)*X,'r');
disp(['g = ',num2str(coef(2)),'  /  sigma^2-g*mu = ',num2str(coef(1))]); 
disp('à comparer aux valeurs attendues: g=0.096 / sigma^2-g*mu=-19.92')


% EOS350 1600ISO

clear all, close all
Ima=double(imread('mire_EOS350D_R1600.png'));
figure, imagesc(Ima), colormap(gray)
moyenne=[];
variance=[];

Row={100:150,200:280,330:380};
Col={70:120,180:240,290:350,410:460,520:580};

for x=1:numel(Row),
    for y=1:numel(Col),
        bloc=Ima(Row{x},Col{y});
        moyenne=[moyenne; mean(bloc(:))];
        variance=[variance;var(bloc(:))];
    end
end

figure, plot(moyenne, variance,'x'), hold on
[coef,coefint]=regress(variance,[ones(length(variance),1),moyenne]);
X=[min(moyenne):10:max(moyenne)];
plot(X,coef(1)+coef(2)*X,'r');
disp(['g = ',num2str(coef(2)),'  /  sigma^2-g*mu = ',num2str(coef(1))]); 
disp('à comparer aux valeurs attendues: g=1.531 / sigma^2-g*mu=-355.3')


% EOS6D 400ISO

clear all, close all
Ima=double(imread('mire_EOS6D_R400.png'));
figure, imagesc(Ima), colormap(gray)
% regardez à la loupe: poussière sur le patch en haut à gauche
moyenne=[];
variance=[];

Row={60:160,220:320,420:520,590:690};
Col={50:150,220:320,410:510,580:680,790:790,940:1040};

for x=1:numel(Row),
    for y=1:numel(Col),
        if (x==1)&(y==numel(Col)) continue; end;
        bloc=Ima(Row{x},Col{y});
        moyenne=[moyenne; mean(bloc(:))];
        variance=[variance;var(bloc(:))];
    end
end

figure, plot(moyenne, variance,'x'), hold on
[coef,coefint]=regress(variance,[ones(length(variance),1),moyenne]);
X=[min(moyenne):10:max(moyenne)];
plot(X,coef(1)+coef(2)*X,'r');
disp(['g = ',num2str(coef(2)),'  /  sigma^2-g*mu = ',num2str(coef(1))]); 
disp('à comparer aux valeurs attendues: g=0.705 / sigma^2-g*mu=-1442')


% EOS6D 6400ISO

clear all, close all
Ima=double(imread('mire_EOS6D_R6400.png'));
figure, imagesc(Ima), colormap(gray)
moyenne=[];
variance=[];

Row={90:190,270:370,450:550,630:730};
Col={70:170,250:350,430:530,620:720,790:890,970:1070};

for x=1:numel(Row),
    for y=1:numel(Col),
        bloc=Ima(Row{x},Col{y});
        moyenne=[moyenne; mean(bloc(:))];
        variance=[variance;var(bloc(:))];
    end
end

figure, plot(moyenne, variance,'x'), hold on
[coef,coefint]=regress(variance,[ones(length(variance),1),moyenne]);
X=[min(moyenne):10:max(moyenne)];
plot(X,coef(1)+coef(2)*X,'r');
disp(['g = ',num2str(coef(2)),'  /  sigma^2-g*mu = ',num2str(coef(1))]); 
disp('à comparer aux valeurs attendues: g=10.52  /  sigma^2-g*mu=-21131')




% Partie II

clear all, close all;

% expérience 1
M=0; V=0; niter=100000;
for i=1:niter; 
    C=dct2(2+3*randn(8,8)); 
    M=M+C/niter;
    V=V+C.^2/niter; 
end; 
M
V-M.^2


% expérience 2
Ima=imread('computernonoise.png');
figure, imshow(Ima);
Bloc=Ima(361:392,201:232);
figure, imshow(Bloc);

% indicatrice HF
T=40;
tima=32;
Q=zeros(tima,tima);
for i=1:tima,
    for j=1:tima,
        if i+j>T,
    	    Q(i,j)=1;
        end;
    end;
end;
Q=boolean(Q);
figure, imagesc(Q);

Bloc=double(Bloc);
sigma=5;
Blocnoise=Bloc+sigma*randn(size(Bloc));
dctblocnoise=dct2(Blocnoise);
disp(['variance estimée: ',num2str(var(dctblocnoise(Q))),' à comparer à ',num2str(sigma^2)])
%mean(dctblocnoise(not(Q)).^2)



% estimation du bruit hétéroscédastique sur une image

clear all, close all

% image synthétique:
Ima=double(imread('computernonoise.png'));
% exemple (à changer)
g=2;
moy=40; 
sigma=5;
Imanoise=g*poissrnd(Ima/g)+moy+sigma*randn(size(Ima));
% equation attendue: g*y + sigma^2-g*moy 
figure, imshow(uint8(Imanoise))
disp(['valeurs attendues: g=',num2str(g),'  /  sigma^2-g*mu=',num2str(sigma^2-g*moy)])
ColomBuades_estimation(Imanoise)


% alternative: image réelle (choisir une des images)
Imanoise=imread('mire_EOS350D_R100.png');
Imanoise=imread('mire_EOS350D_R1600.png');
Imanoise=imread('mire_EOS6D_R400.png');
Imanoise=imread('mire_EOS6D_R6400.png');
Imanoise=imread('clutter_EOS350D_R100.png');
Imanoise=imread('clutter_EOS350D_R1600.png');
Imanoise=imread('clutter_EOS6D_R400.png');
Imanoise=imread('clutter_EOS6D_R6400.png');

figure, imagesc(Imanoise)
ColomBuades_estimation(Imanoise)


