output=magickernel(image, varargin)

function output=magickernel(image, varargin)
% Magic Kernel enlarging
%   MAGICKERNEL(I) doubles the image I size with magic kernel.
%   The enlarging can be repeated to get higher magnification.
%
%   MAGICKERNEL(I, 'DOWN') halves the image I size with magic kernel.
%   The shrinking can be repeated to get smaller image.
%
%   Reference with method comparison:
%       http://johncostella.webs.com/magic/
%
%   Example:
%       imshow(magickernel(imread('eight.tif')))
%
%   See also IMRESIZE.

%   Contributed by Jan Motl ([e-mail])
%   $Revision: 1.0 $  $Date: 2013/02/16 16:58:01 $
%   Note: The enlarging and shrinking subroutines are independant
%         for easy copy&paste. 

% Check parameters
numvarargs = length(varargin);
if numvarargs > 1
    error('myfuns:somefun2Alt:TooManyInputs', ...
        'requires at most 1 optional input');
end
 
optargs = {'up'};       % set defaults for optional inputs
 
optargs(1:numvarargs) = varargin;
[scale] = optargs{:};    % Place optional args in memorable variable names

% Initialization
[height, width, colors] = size(image);

% Perform enlarging/shrinking color by color
if strcmp(scale, 'up') % enlarging
    output = zeros(2*height, 2*width, colors);

    for color = 1:colors
        output(:,:,color) = enlargeOneColorDimension(image(:,:,color));
    end
else                   % shrinking
    output = zeros(ceil(height/2), ceil(width/2), colors);

    for color = 1:colors
        output(:,:,color) = shrinkOneColorDimension(image(:,:,color));
    end
end

% Make sure the output is of the same variable type as the input
output = cast(output, class(image));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Subroutines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% Enlarge black-and-white image %%%%%%%%%%%%%%%%%%%%%%%
    function output = enlargeOneColorDimension(image)
        % Pad the image to avoid dark borders
        image = padarray(image, [1 1], 'replicate', 'both');
        
        % Interlace the image with black strips, vertically and horizontaly
        % Multiplication of sparse matrices is faster than for cycles
        [height, width] = size(image);
        
        matrix1 = spalloc(2*height, height, height);
        for i=1:height
            matrix1(2*i, i) = 1;    %#ok
        end
        
        matrix2 = spalloc(width, 2*width, width);
        for i=1:width
            matrix2(i, 2*i) = 1;    %#ok
        end
        
        output = matrix1 * double(image) * matrix2;
        
        % Perform convolution with magic kernel h
        h = [1/4, 3/4, 3/4, 1/4];
        output = imfilter(output, h'*h);
        
        % Trim the result
        output = output(3:end-2, 3:end-2);
    end

%%%%%%%%%%%%%%%%%%%% Downsample black-and-white image %%%%%%%%%%%%%%%%%%%%%
    function output = shrinkOneColorDimension(image)
        % Pad the image to avoid dark borders
        image = padarray(image, [1 1], 'replicate', 'both');
                
        % Convolve the image with magic kernel to remove high frequencies
        % that could cause moire patterns
        h = 0.5*[1/4, 3/4, 3/4, 1/4];
        output = imfilter(image, h'*h);
        
        % Downsample the image (remove each second row and column)
        [height, width] = size(image);
        height2 = ceil(height/2);
        width2 = ceil(width/2);
        height2l = floor(height/2);
        width2l = floor(width/2);
        
        matrix1 = spalloc(height2, height, height2);
        for j=1:height2l
            matrix1(j, j*2) = 1;    %#ok
        end
        
        matrix2 = spalloc(width, width2, width2);
        for j=1:width2l
            matrix2(j*2, j) = 1;    %#ok
        end
        
        output = matrix1 * double(output) * matrix2;
        
        % Trim the result
        output = output(1:end-1, 1:end-1);
    end

end