function [  ] = perona_malik_1D(  )
% This function demonstrates how to implement the perona-malik algorithm in
% 1D. Explicit, implicit and semi-implicit solutions are shown.

%set algorithm parameters
max_it = 300;
alpha = 2;       %exponent in conduction coefficient
K = 0.02;        %denominator in conduction coefficient
lambda = 0.75;   %time step

%synthesize function to denoise
n = 512;        %number of samples
u0 = zeros(n,1);
u0(ceil(n/3): ceil(2*n/3)) = 1;

figure;
subplot(2,3,1);
plot(u0);
axis([0, n, -0.5, 1.5]);
title('Ground Truth');

% add noise to the function
u0 = u0 + 0.125*randn(n,1);

subplot(2,3,2);
plot(u0);
axis([0, n, -0.5, 1.5]);
title('Noisy Data');

e = ones(n,1);
% Dfx computes forward differences in the x-direction
Dfx = spdiags([-e e],[0:1], n, n);
% Dbx computes backward differences in the x-direction
Dbx = spdiags([-e e],[-1:0], n, n);

I = speye(n,n);

u_exp = u0;
u_imp = u0;
u_mix = u0;

for i = 1:max_it
    i % print iteration number
    
    % Compute conduction coefficients and place on the main diagonal of
    % a sparse matrix. These matrices can be used to scale a vector of
    % derivatives
    Cfx = spdiags(1 ./ (1 + (abs(Dfx*u_exp) ./K).^alpha), 0, n, n);
    Cbx = spdiags(1 ./ (1 + (abs(Dbx*u_exp) ./K).^alpha), 0, n, n);
    Cux = spdiags(1 ./ (1 + (abs(Dux*u_exp) ./K).^alpha), 0, n, n);
    Cdx = spdiags(1 ./ (1 + (abs(Ddx*u_exp) ./K).^alpha), 0, n, n);
    
    % iterate using explicit formulation
    %
    %   u(t+1) = u(t) + lambda*(Cfx*Dfx*u(t) - Cbx*Dbx*u(t))
    %
    A = I + lambda*(Cfx*Dfx + Cbx*Dbx + Cux*Dux + Cdx*Ddx);
    u_exp = A*u_exp;
    
    % iterate using implicit formulation (derivatives applied to u(t+1))
    %
    %   u(t+1) = u(t) + lambda*(Cfx*Dfx*u(t+1) - Cbx*Dbx*u(t+1))
    %
    B = I - lambda(Cfx*Dfx + Cbx*Dbx + Cux*Dux + Cdx*Ddx);
    u_imp = B\u_imp;
      
    %plot results
    subplot(2,3,4);
    plot(u_exp);
    axis([0, n, -0.5, 1.5]);
    title('PM explicit');
    
    subplot(2,3,5);
    plot(u_imp);
    axis([0, n, -0.5, 1.5]);
    title('PM implicit');
    
    subplot(2,3,6);
    plot(u_mix);
    axis([0, n, -0.5, 1.5]);
    title('PM mixed');
    drawnow;
end


end