Laplacian of Guassians Edge Detection

Hello, I'm trying to implement a Laplacian of Guassians Edge Detection using the following code, the issue I have is that Indexing errors pop up when I try to implement the following algorithm

What im trying to do is For each pixel location where the convolution is positive:

a. Test to see of the value of convolution immediately to the left is negative. If so then a zero crossing is present. If the sum of the absolute of the two values is greater than t then this pixel is an edge pixel.

b. repeat step (a) for pixels to the right, above and below the current one.

The edge image, the gradient and filter gradient image should be returned in variables E, F and G respectively

My code:

function [E,F,G] = log_edge(I,N)
% inputs : I - image, N - size of filter 
% outputs: E - edge image, F - filter
if (nargin<2)
  N=5;
end

% force I to be a NxN real number array, and create the laplacian of gaussians filter. Note that F sums to zero has has total energy = 1.
I=double(I(:,:,1));
if (N<=3)
  F=[0 1 0; 1 -4 1; 0 1 0]/8;
else
  F = fspecial('log',N,floor((N-1)/3)/2);
end

% TO PREVENT THE THRESHOLD GOING TOO LOW FOR SMALL FILTER SIZES 3,5,7,9 etc 
threshK=max(1,-0.5*N+7.5); % multiply threshold by this factor


% 1. Create and emtpy array E and and an array G containing the filtered image.
E = zeros(size(I));
G = conv2(I,F,'same');

% 2. compute threshold t (0.75*mean(G)) of the LoG image stored in G (and multiply by threshK)
t = 0.75*mean(G(:))*threshK;

% 3. identify the zero crossing points 
[r, c] = size(I);
for i= 1:r
    for j = 1:c
        if I(i,j-1)>0
% 4. preserve those zero crossing points where the sum of the magnitudes of G accross the zero crossing is > t
           res = abs(j+(j-1));
          if res>t
           E(i,j) = res;
          end
       end
       
       if I(i-1,j)>0
           res = abs(j+(j-1));
           if res>t
           E(i,j) = res;               
           end
       end
       
       if I(i,j+1)>0
           res = abs(j-(j-1));
           if res>t
           E(i,j) = res;
           end
       end
       
       if I(i+1,j)>0
           res = abs(j+(j-1));
           if res>t
           E(i,j) = res;           
           end
       end
       
        
    end
end
return

For reference I have a test script that checks my solution, the script is as follows :

I = double(imread('pout.tif'))/255;

% 13x13 filter
[E,B,G] = log_edge(I,13);
figure
subplot(2,2,1);
imagesc(I); title('Image'); axis equal tight; 
subplot(2,2,2);
imagesc(G); title('LoG Output'); axis equal tight; 
subplot(2,2,3);
imagesc(B); axis equal tight; title('LoG Filter');
subplot(2,2,4);
imagesc(E); title('Edges'); axis equal tight; 
colormap(gray);
drawnow;

% 9x9 filter
[E,B,G] = log_edge(I,15);
figure
subplot(2,2,1);
imagesc(I); title('Image'); axis equal tight; 
subplot(2,2,2);
imagesc(G); title('LoG Output'); axis equal tight; 
subplot(2,2,3);
imagesc(B); axis equal tight; title('LoG Filter');
subplot(2,2,4);
imagesc(E); title('Edges'); axis equal tight; 
colormap(gray);
drawnow;

Thank you for the help and have a good day.

Answer 1

[Anusha Mohan]

I = double(imread('pout.tif'))/255;

% 13x13 filter

[E,B,G] = log_edge(I,13);

figure

subplot(2,2,1);

imagesc(I); title('Image'); axis equal tight;

subplot(2,2,2);

imagesc(G); title('LoG Output'); axis equal tight;

subplot(2,2,3);

imagesc(B); axis equal tight; title('LoG Filter');

subplot(2,2,4);

imagesc(E); title('Edges'); axis equal tight;

colormap(gray);

drawnow;