%Based on the paper http://arxiv.org/abs/1111.1567 
%Options are preset = 0 => discrete time smoothlife
%            preset = 1 => continuous time smoothlife
global M N fourier_kernel_m fourier_kernel_n b1 b2 d1 d2 alphan alpham using1 is_continuous
preset = 1;
NX = 2^8;
NY = NX;
td = 10;
dt = 1/td;
T = 1000;

if preset %Option for continuous time smoothlife
ra = 10;
ri = ra/3;
b = 1;
b1 = 0.257;
b2 = 0.336;
d1 = 0.365;
d2 = 0.551;
alphan = 0.028;
alpham = 0.147;
rk = 1;
using1 = 0;
is_continuous = 1;
else %Option for discrete time smoothlife
ra = 10;
ri = ra/3;
b = 1;
b1 = 0.278;
b2 = 0.365;
d1 = 0.267;
d2 = 0.445;
alphan = 0.028;
alpham = 0.147;
rk = 0;
using1 = 1;
is_continuous = 0;
end

%used to generate antialiased disc shape
function val = linear_step(x, b, radius)
    val = (radius + b/2 - x)/b;
    val(x < radius - b/2) = 1;
    val(x > radius + b/2) = 0;
end
%a sigmoid function
function val = sigma1(x,a,alpha)
    val = 1 ./ (1 + exp(-(x-a)*4/alpha));
end
%a smooth step function
function val = sigma2(x,a,b,alpha)
    val = sigma1(x,a,alpha).*(1 - sigma1(x,b,alpha));
end
%used to mix values from x and y
function val = sigmam(x,y,m,alpha)
    s1 = sigma1(m,0.5,alpha);
    val = x.*(1 - s1) + y.*s1;
end
%variant of the time stepping function, used in discrete time version
function val = s_1(n,m,b1,b2,d1,d2,alphan,alpham)
    val = sigma2(n,sigmam(b1,d1,m,alpham),sigmam(b2,d2,m,alpham),alphan);
end
%variant of the time stepping function, used in continuous time version
function val = s_2(n,m,b1,b2,d1,d2,alphan,alpham)
    val = sigmam(sigma2(n,b1,b2,alphan), sigma2(n,d1,d2,alphan),m,alpham);
end
%if the values go off bounds, return them back to bounds
function val = truncate(x)
    val = x;
    val(x < 0) = 0;
    val(x > 1) = 1;
end
%calculated the next frame in the discrete case or the derivative in the continuous case
function val = fdot(f)
    global M N fourier_kernel_m fourier_kernel_n b1 b2 d1 d2 alphan alpham using1 is_continuous
    fourier_f = fft2(f);
    %applying a convolution of f with m and n with FFT to speed it up
    fourier_m = fourier_f .* fourier_kernel_m;
    fourier_n = fourier_f .* fourier_kernel_n;
    %returning to space domain and normalization
    m = real(ifft2(fourier_m)) / M;
    n = real(ifft2(fourier_n)) / N;
    %applying time stepping functions
    if using1
       val = s_1(n,m,b1,b2,d1,d2,alphan,alpham);
    else
       val = s_2(n,m,b1,b2,d1,d2,alphan,alpham);
    end
    %differential version must return values in range [-1,1]
    if is_continuous
       val = 2*val - 1;
    end
end

[XX, YY] = meshgrid(1:NX,1:NY);
R = sqrt((XX - NX/2 - 1).^2 + (YY - NY/2 - 1).^2);
radius_i = ri;
radius_a = ra;
inner_disc = linear_step(R, b, radius_i);
outer_disc = linear_step(R, b, radius_a);
%calculated kernels for self and neighborhood
kernel_m = inner_disc;
kernel_n = outer_disc - inner_disc;
%integrating to get the normalization constants
M = sum(sum(kernel_m));
N = sum(sum(kernel_n));
%kernels converted to frequency domain to make convolution faster
fourier_kernel_m = fft2(fftshift(kernel_m));
fourier_kernel_n = fft2(fftshift(kernel_n));
%I found that the following initial condition works well for both versions
f = rand(NX,NY) .* (XX < NX/2) .* (YY < NY/2);
f(f < 0.5) = 0;
f(f > 0.5) = 1;
for t = 0:dt:T
    fd = fdot(f);
    if is_continuous
       if rk % applying fourth order Runge-Kutta
          k1 = fd;
          k2 = fdot(truncate(f + k1*dt/2));
          k3 = fdot(truncate(f + k2*dt/2));
          k4 = fdot(truncate(f + k3*dt  ));
          f = f + (k1 + 2*k2 + 2*k3 + k4)*dt/6;
       else % applying euler's method
          f = f + dt*fd;
       end
    else % discrete timestepping does not involve a differential equation
        f = fd;
    end
    f = truncate(f); %keeping the values in [0,1]

    if !is_continuous || mod(t*td, td/2) == 0
       imagesc(f, [0 1]);
       colormap(gray); %plotting in elegant grayscale
       title(sprintf('iteration number %g', t));
       drawnow;
    end
end
