function imageOutput = stackBlur(imageFileName, kernelRadii)
% MATLAB implementation of stack blur by cabbage <stackblur@xcel.me>
%   
%   May/2018: code is finished, usable but rather slow.
%   
%   This function can be applied to RGB images only.
%   The colorChannel thing should be removed if you want to use on
%   gray scale images.
%
%   Stack blur is a fast, almost Gaussian blur, but less resource hungry
%   and more CPU friendly. For algorithm details, check out the tutorial
%   below by Quasimondo:
%   github.com/Quasimondo/QuasimondoJS/blob/master/blur/stackblur.pdf
%   
%   This piece of code is of poor quality and is quite ugly 
%   due to my time pressure. Don't forget to email me if 
%   you have some improvements. :)

    imageInput = double(imread(imageFileName));
    imageOutput = zeros(size(imageInput), 'double');
    kernelRadii = double(round(kernelRadii));
    imageSize = size(imageInput);
    if kernelRadii < 2
        error('kernelRadii must be no less than 2');
    else
        if kernelRadii > min([imageSize(1), imageSize(2)])
            error('kernelRadii must be no larger than min dimension of picture.');
        end
    end
    kernelHorizGrid = zeros(1, 3*imageSize(1), 'double');
    kernelVerGrid = zeros(3*imageSize(2), 1, 'double');
    kernel = double(0);
    meanRate = 1/kernelRadii/kernelRadii;
    
    function horizKernelInitializer(row, colorChannel)
        % Start from the first pixel in a row
        kernelHorizGrid = [imageInput(row,:,colorChannel), imageInput(row,:,colorChannel), imageInput(row,:,colorChannel)];
        kernel = double(0);
        for ii = 1:kernelRadii
            kernel = kernel + kernelHorizGrid(imageSize(2)+ii)*(kernelRadii+1-ii);
        end
        for ii = 1:(kernelRadii-1)
            kernel = kernel + kernelHorizGrid(imageSize(2)+1-ii)*(kernelRadii-ii);
        end
    end

    function vertKernelInitializer(col, colorChannel)
        % Start from the first pixel in a col
        kernelVerGrid = [imageOutput(:,col,colorChannel); imageOutput(:,col,colorChannel); imageOutput(:,col,colorChannel)];
        kernel = double(0);
        for ii = 1:kernelRadii
            kernel = kernel + kernelVerGrid(imageSize(1)+ii)*(kernelRadii+1-ii);
        end
        for ii = 1:(kernelRadii-1)
            kernel = kernel + kernelVerGrid(imageSize(1)+1-ii)*(kernelRadii-ii);
        end
    end

    function horizKernelReloader(col)
        % Move 1 pixel forward
        for ii = 1:kernelRadii
            kernel = kernel + kernelHorizGrid(imageSize(2)+col+ii);
            kernel = kernel - kernelHorizGrid(imageSize(2)+col+1-ii);
        end
    end

    function vertKernelReloader(row)
        % Move 1 pixel downward
        for ii = 1:kernelRadii
            kernel = kernel + kernelVerGrid(imageSize(1)+row+ii);
            kernel = kernel - kernelVerGrid(imageSize(1)+row+1-ii);
        end
    end

    function worker(colorChannel)
        for row = 1:imageSize(1)
            % horizontal one dimensional "convolution"
            horizKernelInitializer(row, colorChannel);
            for col = 1:imageSize(2)
                imageOutput(row, col, colorChannel) = kernel*meanRate;
                horizKernelReloader(col);
            end
        end
        
        for col = 1:imageSize(2)
            % vertical one dimensional "convolution"
            % use the image we already get after horizontal "convolution"
            vertKernelInitializer(col, colorChannel);
            for row = 1:imageSize(1)
                imageOutput(row, col, colorChannel) = kernel*meanRate;
                vertKernelReloader(row);
            end
        end
    end




    for color = 1:imageSize(3)
        worker(color);
    end

    imageOutput = uint8(round(imageOutput));
end