clc
clear
close all

tic
disp('starting...')

x = importdata('input.txt');
z = expandZ4(cell2mat(x)=='#');

for i=1:6
    s       = sumActiveNeighbours(z);
    zNew1   =  z & ((s==2)|(s==3));
    zNew2   = ~z & (s==3);
    z       = zNew1 | zNew2;
    z       = expandZ4(z);
end
disp(sum(sum(sum(sum(z)))))

toc
disp('Finished')

%functions
function z = expandZ4(z)
% if any face of the shape contains nonzero elements, this function will
% pad that face with an extra layer of zeros.

if any(any(any(z(:,:,:,1))))
    zNew                     = false(size(z,1),size(z,2),size(z,3),size(z,4)+1);
    zNew(:,:,:,2:size(z,4)+1)= z;
    z                        = zNew;
end

if any(any(any(z(:,:,:,end))))
    zNew                     = false(size(z,1),size(z,2),size(z,3),size(z,4)+1);
    zNew(:,:,:,1:size(z,4))  = z;
    z                        = zNew;
end

if any(any(any(z(1,:,:,:))))
    zNew                     = false(size(z,1)+1,size(z,2),size(z,3),size(z,4));
    zNew(2:size(z,1)+1,:,:,:)= z;
    z                        = zNew;
end

if any(any(any(z(end,:,:,:))))
    zNew                     = false(size(z,1)+1,size(z,2),size(z,3),size(z,4));
    zNew(1:size(z,1),:,:,:)  = z;
    z                        = zNew;
end

if any(any(any(z(:,1,:,:))))
    zNew                     = false(size(z,1),size(z,2)+1,size(z,3),size(z,4));
    zNew(:,2:size(z,2)+1,:,:)= z;
    z                        = zNew;
end

if any(any(any(z(:,end,:,:))))
    zNew                     = false(size(z,1),size(z,2)+1,size(z,3),size(z,4));
    zNew(:,1:size(z,2),:,:)  = z;
    z                        = zNew;
end

if any(any(any(z(:,:,1,:))))
    zNew                     = false(size(z,1),size(z,2),size(z,3)+1,size(z,4));
    zNew(:,:,2:size(z,3)+1,:)= z;
    z                        = zNew;
end

if any(any(any(z(:,:,end,:))))
    zNew                     = false(size(z,1),size(z,2),size(z,3)+1,size(z,4));
    zNew(:,:,1:size(z,3),:)  = z;
    z                        = zNew;
end
end

function s = sumActiveNeighbours(z)
N = getNeighbours4(z);
s = 0;
for i=1:size(N,1)
    s = s + N{i};
end
end

function N = getNeighbours4(z)

N=cell(80,1);
% returns the 80 neighbours(as cubes) in a cell

idx=1;
for dir1=-1:1
    for dir2=-1:1
        for dir3=-1:1
            for dir4=-1:1
                if ~(dir1==0 && dir2==0 && dir3==0 && dir4==0)
                    N{idx} = shiftZ(shiftZ(shiftZ(shiftZ(z,1,dir1),2,dir2),3,dir3),4,dir4);
                    idx=idx+1;
                end
            end
        end
    end
end
end

function zNew = shiftZ(z,dim,dir)
% shift by 1 on dimension dim(1/2/3/4), dir = 1 shifts one way, dir = -1 the other way, dir=0 does nothing;
zNew=z;
switch dim
    case 1
        if dir==1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(1:end-1,:,:,:) = z(2:end,:,:,:);
        elseif dir==-1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(2:end,:,:,:)   = z(1:end-1,:,:,:);
        end
    case 2
        if dir==1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(:,1:end-1,:,:) = z(:,2:end,:,:);
        elseif dir==-1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(:,2:end,:,:)   = z(:,1:end-1,:,:);
        end
    case 3
        if dir==1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(:,:,1:end-1,:) = z(:,:,2:end,:);
        elseif dir==-1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(:,:,2:end,:)   = z(:,:,1:end-1,:);
        end
    case 4
        if dir==1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(:,:,:,1:end-1) = z(:,:,:,2:end);
        elseif dir==-1
            zNew                = false(size(z,1),size(z,2),size(z,3),size(z,4));
            zNew(:,:,:,2:end)   = z(:,:,:,1:end-1);
        end
    otherwise
        error('Incorrect dimension')
end
end