%%
%% Code for demos presented in lecture 4, CNS186 Winter 2000
%%
%% Adapted from:
%% P. Perona and M. Weber - April 1 1999
%% EE/CNS 148 - Spring 1999
%% Edited Jan 16, 2002
%% Edited Jan 17, 2006
%% (c) California Institute of Technology
%% 

%%
%%create 1D image
%%

%% Standard deviation of noise:
NOISE_STD_VEC = [0.2 0.4 0.8 1.6];  %%% <- Change this value in order to test the robustness of the method
SCREEN_SIZE = get(0,'ScreenSize');
WINDOW_POS = round(SCREEN_SIZE*0.95);

for nnn=1:length(NOISE_STD_VEC),
    figure(nnn); clf;
    NOISE_STD = NOISE_STD_VEC(nnn),
    %% A few formatting parameters and initialization of ima array
    LN_WDTH = 2; 	%% Thickness of line in plots
    HALF_IMA_LEN = 60; %% Half length of the 1D image
    IMA_LEN = 2*HALF_IMA_LEN;
    L = [2:IMA_LEN IMA_LEN];
    R = [1 1:(IMA_LEN-1)];
    
    
    %%%%%%%
    %% Generate signal
    %%%%%%%
    ima =  [-1*ones(1,HALF_IMA_LEN) 1*ones(1,HALF_IMA_LEN)]; %% This is a step edge, you may change this
    figure(nnn); subplot(4,1,1); plot(ima,'LineWidth',LN_WDTH);
    title('Image with signal and no noise.');
    axis([1 IMA_LEN -4 4]);
    %axis off;
    
    %% Add noise to image
    ima = ima + NOISE_STD*randn(1,IMA_LEN); %% This is white Gaussian noise
    figure(nnn); set(nnn,'Position',WINDOW_POS,'WindowStyle','modal'); subplot(4,1,2); plot(ima,'LineWidth',LN_WDTH);
    title(['Image with white Gaussian noise, \sigma=' num2str(NOISE_STD)]);
    %axis off;
    
    
    %%%%%%%
    %% Match-filter the image and calculate an appropriate threshold
    %%%%%%%
    
    %% Using a box-filter
    k_box = [zeros(1,8) -1*ones(1,8) 1*ones(1,8) zeros(1,8)];%% Set the kernel to be the matched filter
    %% It needs to be flipped because of how conv2() will use it
    k_box = k_box / norm(k_box); %% Normalize the kernek k
    k_norm = norm(k_box);
    
    out = conv2(ima,-k_box,'same');
    edges = find( (out>out(R))&(out>out(L))); %% Find position of local maxima
    figure(nnn); subplot(4,1,3); plot(out,'LineWidth',LN_WDTH);
    THRESH = 3*NOISE_STD*(k_norm - 0);
    hold on; plot([1 IMA_LEN],[THRESH THRESH],'g','LineWidth',LN_WDTH); hold off;
    hold on; plot([1:length(k_box)],THRESH+(1+k_box)*5,'r','LineWidth',LN_WDTH); hold off;
    hold on; plot(edges,out(edges),'ro','MarkerSize',10); hold off;
    i_detected = find(out(edges)>THRESH);
    hold on; plot(edges(i_detected),out(edges(i_detected)),'r.','MarkerSize',20); hold off;
    %axis off
    title('Output of box filter; the thresh is set at half-point between the E d and E d_n');
    
    
    %% Using a Gaussian filter
    k_gauss = [1 2 1]/4;
    k_smooth = k_box;
    for i=1:8,
        k_smooth = conv2(k_smooth,k_gauss,'same');
    end;
    k_smooth = k_smooth / norm(k_smooth);
    
    out = conv2(ima,-k_smooth,'same');
    edges = find( (out>out(R))&(out>out(L))); %% Find position of local maxima
    figure(nnn); subplot(4,1,4); plot(out,'LineWidth',LN_WDTH);
    THRESH = 3*NOISE_STD*(k_norm - 0);
    hold on; plot([1 IMA_LEN],[THRESH THRESH],'g','LineWidth',LN_WDTH); hold off;
    hold on; plot([1:length(k_smooth)],THRESH+(1+k_smooth)*5,'r','LineWidth',LN_WDTH); hold off;
    hold on; plot(edges,out(edges),'ro','MarkerSize',10); hold off;
    i_detected = find(out(edges)>THRESH);
    hold on; plot(edges(i_detected),out(edges(i_detected)),'r.','MarkerSize',20); hold off;
    %axis off
    title('Output of smooth filter; the thresh is set at half-point between the E d and E d_n');
    
    eval(['print -dpsc edge1D_' num2str(nnn)]);
end