%%
close all
% Read image
image = imread('toronto_original.png');
[rows, columns, channels] = size(image);

figure('Name', 'Original');
imshow(image)

for col = 1:columns
    for row = 1:rows
        if mod(row, 2) && mod(col, 2) %Both odd
            % Blue
            image(row,col,[1,2]) = 0;
        elseif ~mod(row, 2) && ~mod(col, 2) %Both even
            % Red
            image(row,col,[2,3]) = 0;
        else %One odd one even
            % Green
            image(row,col,[1,3]) = 0;
        end
    end
end

figure('Name', 'CFA')
imshow(image)

%%

BIimage = zeros(rows, columns, channels, 'uint8');

for col = 1:columns
    for row = 1:rows
        % Accounts for image edges
        rowMinus = max(row-1, 1);
        rowPlus = min(row+1, rows);
        colMinus = max(col-1, 1);
        colPlus = min(col+1, columns);

        submatrix = image(rowMinus:rowPlus,colMinus:colPlus,:);

        if mod(row, 2) && mod(col, 2) %Both odd
            % Blue
            red = mean( nonzeros( submatrix(:,:,1) ));
            green = mean( nonzeros( submatrix(:,:,2) ));
            blue = image(row,col,3);

        elseif ~mod(row, 2) && ~mod(col, 2) %Both even
            % Red
            blue = mean( nonzeros( submatrix(:,:,3) ));
            green = mean( nonzeros( submatrix(:,:,2) ));
            red = image(row,col,1);

        else %One odd one even
            % Green
            red = mean( nonzeros( submatrix(:,:,1) ));
            blue = mean( nonzeros( submatrix(:,:,3) ));
            green = image(row,col,2);
        end
        BIimage(row,col,1) = red;
        BIimage(row,col,2) = green;
        BIimage(row,col,3) = blue;
    end
end

figure('Name', 'Bilinear Interpolation')
imshow(BIimage)