type BiSmooth.m

% smooth image A using a binomial filters of order p (should be even)
function [SI] = BiSmooth(A,p);

X = [1 2 1];
c = 4;

for i=2:2:p-2,
	  X = conv(X,[1 2 1]);
  c = c*4;
end;

X, c

I1 = filter2(X/c,A);
I2 = filter2(X'/c,I1);
SI = uint8(I2);


I = imread('background.tif');


figure, subplot(2,2,1); imshow(I);
Ig = rgb2gray(I); 
subplot(2,2,1); imshow(Ig);

Ig2 = BiSmooth(Ig,2);

X =

     1     2     1


c =

     4

subplot(2,2,2); imshow(Ig2);
Ig4 = BiSmooth(Ig,4);

X =

     1     4     6     4     1


c =

    16

subplot(2,2,3); imshow(Ig4);
Ig16 = BiSmooth(Ig,16);

X =

  Columns 1 through 8

           1          16         120         560        1820        4368        8008       11440

  Columns 9 through 16

       12870       11440        8008        4368        1820         560         120          16

  Column 17

           1


c =

       65536

subplot(2,2,4); imshow(Ig16);
figure, plot(Ig(240,:));
hold;
Current plot held
plot(Ig2(240,:),'r');
plot(Ig4(240,:),'g');
plot(Ig16(240,:),'c');


type Gradient.m

% Partial derivatives in X and Y directions
function [Gx,Gy] = Gradient(A, p);

X = [-1 0 1]/2; % filter
Y = [1 0 -1]'/2; % filter

Gx = filter2(X,A); % correlation
Gy = filter2(Y,A);
[M,N] = size(A);
Gx(1:p+1,:) = 0;
Gx(M-p-1:M,:) = 0;
Gx(:,1:p+1) = 0;
Gx(:,N-p-1:N) = 0;
Gy(1:p+1,:) = 0;
Gy(M-p-1:M,:) = 0;
Gy(:,1:p+1) = 0;
Gy(:,N-p-1:N) = 0;





[Gx,Gy] = Gradient(Ig,0);



m = sqrt(Gx.^2 + Gy.^2);



max(max(m))

ans =

  127.8173



figure, imshow(uint8(m*2)); 
[Gx2,Gy2] = Gradient(Ig,2);
m2 = sqrt(Gx2.^2 + Gy2.^2);
max(max(m2))

ans =

  108.9002

m2 = max(max(m2))

m2 =

  108.9002

figure, imshow(uint8(255*m/m2));
title('p = 2'); 
[Gx4,Gy4] = Gradient(Ig,4);
m4 = sqrt(Gx4.^2 + Gy4.^2);
m2 = sqrt(Gx2.^2 + Gy2.^2);

figure, imshow(uint8(255*m2/(max(max(m2)))));
title('p = 2');
[Gx4,Gy4] = Gradient(Ig,4);
m4 = sqrt(Gx4.^2 + Gy4.^2);
figure, imshow(uint8(255*m4/(max(max(m4)))));
type Gradient.m

% Partial derivatives in X and Y directions
function [Gx,Gy] = Gradient(A, p);

X = [-1 0 1]/2; % filter
Y = [1 0 -1]'/2; % filter

Gx = filter2(X,A); % correlation
Gy = filter2(Y,A);
[M,N] = size(A);
Gx(1:p+1,:) = 0;
Gx(M-p-1:M,:) = 0;
Gx(:,1:p+1) = 0;
Gx(:,N-p-1:N) = 0;
Gy(1:p+1,:) = 0;
Gy(M-p-1:M,:) = 0;
Gy(:,1:p+1) = 0;
Gy(:,N-p-1:N) = 0;





[Gx2,Gy2] = Gradient(Ig2,2);
[Gx4,Gy4] = Gradient(Ig4,4);
[Gx16,Gy16] = Gradient(Ig16,16);
m2 = sqrt(Gx2.^2 + Gy2.^2);
m4 = sqrt(Gx4.^2 + Gy4.^2);
m16 = sqrt(Gx16.^2 + Gy16.^2);
figure, imshow(uint8(255*m2/(max(max(m2)))));
figure, imshow(uint8(255*m4/(max(max(m4)))));
figure, imshow(uint8(255*m16/(max(max(m16)))));
[Gx,Gy] = Gradient(Ig,0);
m = sqrt(Gx.^2 + Gy.^2);


figure, plot(m(240,:)); hold;
Current plot held
plot(m2(240,:),'r'); plot(m4(240,:),'g'); plot(m16(240,:),'c'); 



type NMS2.m

% This is vectorized version

% Non-maximum suppression for a gradient image;
% magnitudes under t are ignored
% returns a binary image where maximal gradient values are non-zero
function [m1] = NMS2(Gx,Gy,m,t);

[M,N] = size(Gx);

% m = sqrt(Gx.^2+Gy.^2);  % magnitude
I1 = find(m<t);
m(I1) = 0;    % zero values under t

[I,J] = find(m>0);
aGx = abs(Gx);
aGy = abs(Gy);

Z = zeros([M,N]);  % binary image

for i=1:size(I),
	if I(i)>1 & I(i)<M & J(i)>1 & J(i)<N
           if Gx(I(i),J(i))>0 
              dx = 1;
           else 
              dx = -1;
           end;
           if Gy(I(i),J(i))>0 
              dy = 1;
           else
	      dy = -1;
           end;

        if aGx(I(i),J(i))>aGy(I(i),J(i)) % x is dominant direction
	   m01 = m(I(i),J(i)+dx);
           m11 = m(I(i)-dy,J(i)+dx);
           m_01 = m(I(i),J(i)-dx);
           m_11 = m(I(i)+dy,J(i)-dx);
           m1 = (aGx(I(i),J(i))-aGy(I(i),J(i)))*m01 + aGy(I(i),J(i))*m11;
           m2 = (aGx(I(i),J(i))-aGy(I(i),J(i)))*m_01 + aGy(I(i),J(i))*m_11;
           mag = m(I(i),J(i))*aGx(I(i),J(i));
           if (mag>m2) & (mag>=m1) 
              Z(I(i),J(i)) = 1;
           end
	else % y is dominant direction
	   m01 = m(I(i)-dy,J(i));
           m11 = m(I(i)-dy,J(i)+dx);
           m_01 = m(I(i)+dy,J(i));
           m_11 = m(I(i)+dy,J(i)-dx);
           m1 = (aGy(I(i),J(i))-aGx(I(i),J(i)))*m01 + aGx(I(i),J(i))*m11;
           m2 = (aGy(I(i),J(i))-aGx(I(i),J(i)))*m_01 + aGx(I(i),J(i))*m_11;
           mag = m(I(i),J(i))*aGy(I(i),J(i));
           if (mag>m2) & (mag>=m1) 
              Z(I(i),J(i)) = 1;
	   end
	end
      end % if
end % for

I1 = find(Z==0);  % only keep points where Z=1
m(I1)=0;
m1 = uint8(m);


type Hysteresis.m

% Hysteresis thresholding: keep all points with m>t connected to 
% points with m>t2 (=2.5t)
function [m1] = Hysteresis(m, t)

[L, numc] = bwlabel(m);  
t2 = 2.5*t;

[I,J] = find((L>0) & (m>t2));  % find labeled points with mag>t2
comps = zeros(numc,1);  % initialize comps

for i=1:size(I),        % components with at least one point mag.>t2
  comps(L(I(i),J(i))) = 1;
end

[I,J] = find(L>0);         % find non-zero points
for i=1:size(I),       
  L(I(i),J(i)) = comps(L(I(i),J(i)));  %replace connected components labels
end;
   
m1 = zeros(size(m));     % initialize m
I2 = find(L>0);          % L can be treated as a vector
m1(I2) = m(I2);

% for i=1:numc,    % remove connected components that are too small
%     i, comps(i)
%     if comps(i)>0
%         ind = find(L == i);
%         i, length(ind)
%         if (length(ind) < MinEdgeLen) 
%             comps(i) = 0;
%         end;
%     end;
% end;
        




type NMS2.m

% This is vectorized version

% Non-maximum suppression for a gradient image;
% magnitudes under t are ignored
% returns a binary image where maximal gradient values are non-zero
function [m1] = NMS2(Gx,Gy,m,t);

[M,N] = size(Gx);

% m = sqrt(Gx.^2+Gy.^2);  % magnitude
I1 = find(m<t);
m(I1) = 0;    % zero values under t

[I,J] = find(m>0);
aGx = abs(Gx);
aGy = abs(Gy);

Z = zeros([M,N]);  % binary image

for i=1:size(I),
	if I(i)>1 & I(i)<M & J(i)>1 & J(i)<N
           if Gx(I(i),J(i))>0 
              dx = 1;
           else 
              dx = -1;
           end;
           if Gy(I(i),J(i))>0 
              dy = 1;
           else
	      dy = -1;
           end;

        if aGx(I(i),J(i))>aGy(I(i),J(i)) % x is dominant direction
	   m01 = m(I(i),J(i)+dx);
           m11 = m(I(i)-dy,J(i)+dx);
           m_01 = m(I(i),J(i)-dx);
           m_11 = m(I(i)+dy,J(i)-dx);
           m1 = (aGx(I(i),J(i))-aGy(I(i),J(i)))*m01 + aGy(I(i),J(i))*m11;
           m2 = (aGx(I(i),J(i))-aGy(I(i),J(i)))*m_01 + aGy(I(i),J(i))*m_11;
           mag = m(I(i),J(i))*aGx(I(i),J(i));
           if (mag>m2) & (mag>=m1) 
              Z(I(i),J(i)) = 1;
           end
	else % y is dominant direction
	   m01 = m(I(i)-dy,J(i));
           m11 = m(I(i)-dy,J(i)+dx);
           m_01 = m(I(i)+dy,J(i));
           m_11 = m(I(i)+dy,J(i)-dx);
           m1 = (aGy(I(i),J(i))-aGx(I(i),J(i)))*m01 + aGx(I(i),J(i))*m11;
           m2 = (aGy(I(i),J(i))-aGx(I(i),J(i)))*m_01 + aGx(I(i),J(i))*m_11;
           mag = m(I(i),J(i))*aGy(I(i),J(i));
           if (mag>m2) & (mag>=m1) 
              Z(I(i),J(i)) = 1;
	   end
	end
      end % if
end % for

I1 = find(Z==0);  % only keep points where Z=1
m(I1)=0;
m1 = uint8(m);

m1 = NMS2(Gx16,Gy16,m16,2);


figure, imshow(m1);
figure, imshow((uint8(m1/max(max(m1)))));
figure, imshow((uint8(255*m1/max(max(m1)))));
figure, imshow(uint8(255*(m1/max(max(m1))))));
??? figure, imshow(uint8(255*(m1/max(max(m1))))));
                                                |
Error: Unbalanced or unexpected parenthesis or bracket.

figure, imshow(uint8(255*(m1/max(max(m1)))));
m1 = NMS2(Gx16,Gy16,m16,0.5);
figure, imshow(uint8(255*(m1/max(max(m1)))));


figure, imshow(m1);
diary