%% Skript pro vytvoření šumových snímků pro ověřené šumového modelu. 

% Nasimuluje se šum v průměrných snímcích pro každou expozici (ty
% předpokládame, že mají minimum šumu), vezme se z nich obdélník ve středu
% šumového snímku a umístí se do výsledného obrázku. Tohle je zopakováno
% pro každou expozici a šum je uložen od nejmenší expozice po největší. To 
% samé se provede pro šum z naměřených bright framů.

f=@(x) 3.17*10^(-5)*(x)-1.85*10^(-6);
final=zeros(1250,1648);

for i=1:9
    ProjectionData=Data_loading; 
    data=ProjectionData.Data;
    fprintf('Nacten %f snimek. \n', i)
    [data, l]=remove_hot_pixels(data,7);
    data=data/(2^16);
    velikost=size(data);
    pocet=velikost(3);
    a=velikost(1);
    b=velikost(2);
    
    for z=1:pocet % normalizace šumu, jako normální bereme první, postupujeme přes průměr
        % vyplníme outliery tak, aby nám nezasahovaly do průměrného snímku
        filled_data(:,:,z)=filloutliers(data(:,:,z),'center','percentiles',[0.5 99.5]);
        prum(z) = mean(filled_data(:, :, z),'all');                           
        Norm_data(:, :, z) = data(:, :, z) *(prum(1) / prum(z)); % znormalizované snímky
    end
    
    prumer=0; % spočtení průměrného snímku
    for j=1:pocet
        prumer=prumer+Norm_data(:,:,j);
    end
    prumer=prumer/pocet;
    str(i)=mean(prumer(:));
    local_var=f(prumer);
    p=round(b/9);
    A=imnoise(prumer,'localvar',im2double(local_var));
    B=A-prumer;
    variance_simulated(i)=var(B(:));
    
    data=filloutliers(data,"nearest","percentiles",[1.5 98.5]);
    [noise]=Noise_extraction(data);
    
    M(i)=mean(noise(:));
    V(i)=var(noise(:));
    
    
    final_sim(:,(1+p*(i-1)):(i*p))=A(:,p*4:5*p-1);
    
    final(:,(1+p*(i-1)):(i*p))=prumer(:,p*4:5*p-1);
    
    for y=1:10
       X=Norm_data(:,:,y);
       real=X-prumer;
       real=real(:,p*4:5*p-1);
       var_real(i,y)=var(real(:));
    end
    
    
    noise_data(:,(1+p*(i-1)):(i*p))=data(:,p*4:5*p-1,1);

end

if size(final_sim,2)~=size(final,2)

    final_sim(:,end+1)=0;
end

if size(noise_data,2)~=size(final,2)
   noise_data(:,end+1)=0;
end

figure,imshow(final,[])

figure, imshow(noise_data,[])

real_sum=noise_data-final;
for z=1:9
    cast=real_sum(:,(1+p*(z-1)):(z*p));
    var_R(z)=var(cast(:));
end

figure,imshow(final_sim,[])
sim_sum=final_sim-final;
for z=1:9
    cast=sim_sum(:,(1+p*(z-1)):(z*p));
    var_sim(z)=var(cast(:));
end

real_sum=noise_data-final;
figure, imshow(sim_sum,[min(sim_sum(:)), max(sim_sum(:))])
figure, imshow(real_sum,[min(sim_sum(:)), max(sim_sum(:))])
x=[1 5 10 15 20 25 30 35 40];
figure
plot(x,var_sim,'o')


figure
plot(x,var_real(:,1),'.')
hold on 
plot(x,var_real(:,2),'.')
plot(x,var_real(:,3),'.')
plot(x,var_real(:,4),'.')
plot(x,var_real(:,5),'.')
plot(x,var_real(:,6),'.')
plot(x,var_real(:,7),'.')
plot(x,var_real(:,8),'.')
plot(x,var_real(:,9),'.')
plot(x,var_real(:,10),'.')
plot(x, var_sim, 'o')
legend('simulated', 'real');

% imwrite(final,'final.tif')
% imwrite(noise_data, 'real.tif')
% imwrite(A, 'simulated.tif')