function output = NESIF( stack , RS, nbits)
% code for the paper "Measuring the noise of digital imaging sensors 
% by stacking raw images affected by vibrations and illumination flickering"
% SIAM Journal on Imaging Sciences
% vol. 8, no. 1, p. 611-643
% Frederic Sur, Michel Grediac
% 2015
% usage:
% stack: input image stack (size1 x size 2 x nr_im double matrix)
% RS: rolling shutter effect compensation (1=yes, 0=no)
% nbits: number of bits per channel (to detect clipped pixels)
% output: [gain (g), offset (mu), readout / quantization noise (sigma),
% sigma^2-g.mu, flickering st.d., x-vibration st.d., y-vibration st.d.,
% covariance of vibration components]
% Remark: to use with a real image stack, use dcraw to linearly convert the RAW
% files and extract a channel (see paper)

disp(' ')
disp('NESIF v. 1.0')
disp('Please refer to:')
disp('Frédéric Sur and Michel Grédiac')
disp('Measuring the noise of digital imaging sensors by stacking raw images affected by vibrations and illumination flickering')
disp('SIAM Journal on Imaging Sciences')
disp('vol. 8, no. 1, p. 611-643')
disp('2015.')
disp(' ')
drawnow('update')

nr_im=size(stack,3);
image_moyenne=mean(double(stack),3);
var_bruit=var(double(stack),0,3);
seuil=2^nbits-1;

inliers_1=find(all(stack<seuil,3));
[GX,GY]=gradient(image_moyenne);
[var_bruit_sum_x,var_bruit_sum_y]=var_grad_comp(var_bruit);    
inliers_2=find(chi2cdf(nr_im*(GX(:).^2./var_bruit_sum_x(:)+GY(:).^2./var_bruit_sum_y(:)),2)<.99);
inliers=intersect(inliers_1,inliers_2); % inliers are not saturated and not along edges

Observ=[image_moyenne(inliers_1), image_moyenne(inliers_1).^2, GX(inliers_1).^2, GY(inliers_1).^2, 2*GX(inliers_1).*GY(inliers_1)];
[b,stats]=robustfit(Observ,var_bruit(inliers_1));
% from eq. (2.20)
sample_stda=sqrt(b(4)/(1+b(3)));
sample_stdb=sqrt(b(5)/(1+b(3)));
sample_covab=b(6)/(1+b(3));

figure, plot(image_moyenne(:),var_bruit(:),'xr'), xlabel('gray level'), ylabel('estimated variance')
title('Variance against average gray level')
hold on
plot(image_moyenne(inliers_1),var_bruit(inliers_1),'xg','LineWidth',2)
plot(image_moyenne(inliers),var_bruit(inliers),'xb','LineWidth',2)
grid on
V=axis;
XX=[V(1):1:V(2)];
YY=b(1)+b(2)*XX+b(3)*XX.^2;
plot(XX,YY,'r--','LineWidth',2)
drawnow;


p=.99;
vpowstar=zeros(1,nr_im);
q99=quantile(image_moyenne(inliers),p);
I=find(image_moyenne(inliers)>=q99);
var_q=min(var_bruit(inliers(I)));
I=find((image_moyenne(inliers)>=q99-3*sqrt(var_q/nr_im))&(image_moyenne(inliers)<=q99+3*sqrt(var_q/nr_im)));  % eq. (3.7)
W=exp(-(image_moyenne(inliers(I))-q99).^2/(2*var_q/nr_im));
W=W/sum(W);
for i=1:nr_im
    image_temp=stack(:,:,i);
    vpowstar(i)=sum(W.*image_temp(inliers(I))./image_moyenne(inliers(I)));  % eq. (3.8)
end;

unplusgammat=zeros(1,1,length(vpowstar));
unplusgammat(:)=vpowstar(:);

stack_b=double(stack)-repmat(unplusgammat,[size(stack,1),size(stack,2),1]).*repmat(image_moyenne,[1,1,nr_im]);

vindstar=zeros(1,size(stack_b,3));
p=0.01;
q01=quantile(image_moyenne(inliers),p);
I=find(image_moyenne(inliers)<=q01);
var_q=max(var_bruit(inliers(I)));
I=find((image_moyenne(inliers)<=q01+3*sqrt(var_q/nr_im))&(image_moyenne(inliers)>=q01-3*sqrt(var_q/nr_im))); % eq (3.13)
W=exp(-(image_moyenne(inliers(I))-q01).^2/(2*var_q/nr_im));
W=W/sum(W);
for i=1:size(stack_b,3)
    image_temp=stack_b(:,:,i);
    vindstar(i)=sum(W.*image_temp(inliers(I)));  % eq. (3.14)
end;
[coef,stats]=robustfit(unplusgammat(:)-1,vindstar(:),'bisquare',4.685,'off');

mud=coef*q99/(coef-q99+q01);  % eq. (3.16)


if (RS==0)
    % no rolling shutter effect
    unplusgammat=1+(unplusgammat-1)*(1+mud/(q99-mud));  % eq. (3.17)
    vargamma=var(unplusgammat);
    stack_b=double(stack)-repmat(unplusgammat,[size(stack,1),size(stack,2),1]).*repmat(image_moyenne-mud,[1,1,nr_im]);  % eq. (3.18)

else
    % if rolling shutter effect (section 3.5)
    % estimation of gamma(t) on a horizontal strip    
    sigma_x=50;
    p=.9;
    vpowstar=zeros(size(stack,1),nr_im);
    q99_tab=zeros(size(stack,1),1);
    I99_tab=zeros(size(stack,1),1);
    [inliersx,inliersy]=ind2sub(size(image_moyenne),inliers);
    for x=1:size(stack,1),
        inliers_c=inliers((inliersx>x-round(sigma_x))&(inliersx<x+round(sigma_x)));
        q99=quantile(image_moyenne(inliers_c),p);
        q99_tab(x)=q99;
        I=find(image_moyenne(inliers)>=q99);
        var_q=min(var_bruit(inliers(I)));
        I=find((image_moyenne(inliers)>=q99-3*sqrt(var_q/nr_im))&(image_moyenne(inliers)<=q99+3*sqrt(var_q/nr_im)));
        [Ix,Iy]=ind2sub(size(image_moyenne),inliers(I));
        W=exp(-(image_moyenne(inliers(I))-q99).^2/(2*var_q/nr_im));
        W2=exp(-(Ix-x).^2/(2*sigma_x^2));
        I99_tab(x)=sum(W2.^2);
        W2=W2/sum(W.*W2);
        for i=1:nr_im
            image_temp=stack(:,:,i);
            vpowstar(x,i)=sum(W2.*W.*image_temp(inliers(I))./image_moyenne(inliers(I)));
        end;
    end;
    
    unplusgammat=zeros(size(stack,1),1,size(stack,3));
    unplusgammat(:)=vpowstar(:);
    for x=1:size(stack,1)
        unplusgammat(x,:,:)=1+(unplusgammat(x,:,:)-1)*(1+mud/(q99_tab(x)-mud));
    end;
    vargamma=var(unplusgammat(:));    
    stack_b=double(stack)-repmat(unplusgammat,[1,size(stack,2),1]).*repmat(image_moyenne-mud,[1,1,nr_im]); 

end;

var_bruit_b=var(double(stack_b),0,3);
Observ=[ones(length(inliers),1), image_moyenne(inliers)];
[b,stats]=regress(var_bruit_b(inliers), Observ);  % from eq. (3.24)
% or:
% [b,stats]=robustfit(image_moyenne(inliers),var_bruit_b(inliers));

output=[b(2), mud, sqrt(b(1)+b(2)*mud), b(1), sqrt(vargamma), sample_stda, sample_stdb, sample_covab];

figure, plot(image_moyenne(:),var_bruit_b(:),'xr'), xlabel('gray level'), ylabel('variance')
title('Corrected variance against average gray level')
hold on
plot(image_moyenne(inliers_1),var_bruit_b(inliers_1),'xg')
plot(image_moyenne(inliers),var_bruit_b(inliers),'xb')
grid on
V=axis;
XX=[V(1):1:V(2)];
YY=b(1)+b(2)*XX;
plot(XX,YY,'r--','LineWidth',2)
drawnow;


end

