function ColomBuades_estimation(Imanoise)
% il s'agit d'une version simplifiée de:
% Miguel Colom, and Antoni Buades,
% Analysis and Extension of the Ponomarenko et al. Method,
% Estimating a Noise Curve from a Single Image,
% Image Processing On Line, 3 (2013), pp. 173–197.
% https://doi.org/10.5201/ipol.2013.45    
  
% indicatrice HF
T=12;
tbloc=8;
Q=zeros(tbloc,tbloc);
for i=1:tbloc,
    for j=1:tbloc,
        if i+j>T,
    	    Q(i,j)=1;
        end;
    end;
end;
Q=boolean(Q);

% indicatrice BF (enlever coef (1,1) qui est toujours la somme des niveaux)
QBF=not(Q);
QBF(1,1)=0;

nblocs=floor(size(Imanoise)/tbloc);
moyenne=[];
varblocHF=[];
varblocBF=[];

% on calcule variance dans partie HF et partie BF
for i=0:nblocs(1)-1,
    for j=0:nblocs(2)-1;
      bloc=Imanoise(tbloc*i+1:tbloc*(i+1),tbloc*j+1:tbloc*(j+1));
      m=mean(bloc(:));  
      moyenne=[moyenne;m];
      dctbloc=dct2(bloc); 
      varblocHF=[varblocHF;var(dctbloc(Q))];
      varblocBF=[varblocBF;var(dctbloc(QBF))]; 
   end;
end;

% régression sur les variances des parties HF, censées correspondre au
% bruit
[coef,coefint]=regress(varblocHF,[ones(length(varblocHF),1),moyenne]);
figure, plot(moyenne,varblocHF,'x');
hold on
X=[min(moyenne):10:max(moyenne)];
plot(X,coef(1)+coef(2)*X,'r')
axis([min(moyenne),max(moyenne),min(varblocHF),quantile(varblocHF,.95)])

disp('estimation 1')
disp(['g = ',num2str(coef(2)),'  /  sigma^2-g*mu = ',num2str(coef(1))]); 


% parmi des blocs dans la même moyenne, on ne garde que les blocs 
% correspondant au 5% de variance BF les plus faibles 
% (ce sont des blocs homogènes, dont la variance HF n'est pas perturbée par la texture)
prop=0.05;
Indices=[];
for val=min(moyenne)+10:20:max(moyenne),
    I=find((moyenne<val)&(moyenne>=val-10));
    [sortvarblocBF,Index]=sort(varblocBF(I));
    Indices=[Indices;I(Index(1:round(numel(I)*prop)))];
end;
moyennebis=moyenne(Indices);
variancebis=varblocHF(Indices);

[coef,coefint]=regress(variancebis,[ones(length(variancebis),1),moyennebis]);
figure, plot(moyenne,varblocHF,'oc',moyennebis,variancebis,'x');
hold on
X=[min(moyenne):10:max(moyennebis)];
plot(X,coef(1)+coef(2)*X,'r')
axis([min(moyenne),max(moyenne),min(varblocHF),quantile(varblocHF,.95)])

disp('estimation 2')
disp(['g = ',num2str(coef(2)),'  /  sigma^2-g*mu = ',num2str(coef(1))]); 


