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);
        
        % 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);
        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);
        u_imp = B\u_imp;
        
        % iterate using mixed explicit/implicit formulation: the average of
        % explicit and implicit formulation
        %
        %          u(t) = A*u(t-1)
        %  +     B*u(t) = u(t-1)
        %  = (B+I)*u(t) = (A+I)*u(t-1)
        %
        u_mix = (B+I)\((A+I)*u_mix);
        
        %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