function [sigma] = noise_estimation(data)
% Noise estimation function based on modified algorith described in:
% J. Immerkaer, “Fast Noise Variance Estimation”, Computer Vision and Image Understanding, Vol. 64, No. 2, pp. 300-302, Sep. 1996
% https://doi.org/10.1006/cviu.1996.0060
% Output: Standard deviation of noise
% modified by Jakub Šalplachta

if size(data, 1)==1 || size(data, 2)==1
    data = squeeze(data);
end

sigma_pom = zeros(1,size(data,3));
mask = zeros(size(data,1),size(data,2));
img = data(:,:,1);
mask(img>0) = 1;
mask = mask(2:end-1,2:end-1);
kernel = gpuArray([1 -2 1; -2 4 -2; 1 -2 1]);
for i=1:size(data,3)
    img = data(:,:,i);
    img = gpuArray(img);
    img_conv = conv2(img, kernel, 'valid');
    img_std = stdfilt(img);
    img_std = img_std(2:end-1,2:end-1);
    img_std = proc.normalize_quick(img_std);
    temp = zeros(size(img_std));
    temp(img_std<(2*mean(img_std(:))+1*std(img_std(:)))) = 1;
    temp = temp.*mask;
    val = (sqrt(0.5*pi)/(6*sum(temp(:))))*sum(sum((abs(img_conv.*temp))));
    sigma_pom(i) = gather(val);
end
sigma = mean(sigma_pom);
end