Homework 5: Matlab Code

Sobel operator
Edge map for Sobel
Laplacian
Double edge map for laplacian
Spatial Average Low-pass filter

spatavg.m

function imOut = spatavg(im, N)
% B = SPATAVG(A,N)		NxN smoothing spatial average filter.

% filter coefficients
A = ones(N,N);
A = A .* (1/N);

[R,C] = size(im);

% pad near the edges with the symmetric extension
if mod(N,2)
   padTL = (N-1)/2;	% left and top
   padBR = (N-1)/2;	% bottom and right
else
   padTL = N/2 - 1;
   padBR = N/2;
end
offset = padTL;
paddedIm = zeros(R+N-1, C+N-1);
[padR,padC] = size(paddedIm);
% copy the main image
paddedImg((offset+1):(offset+R), (offset+1):(offset+C)) = im;
% now copy up into the symmetric extension - top, left, bottom, right
% what a pain - is there a simpler way? (yes-sympad. thanks Lexing. but since
%   this works already, it's easier to leave it... :)
paddedImg(1:padTL,(offset+1):(offset+C)) = flipud(im(1:padTL,:));
paddedImg((offset+1):(offset+R),1:padTL) = fliplr(im(:,1:padTL));
paddedImg((padR-padBR+1):padR,(offset+1):(offset+C)) = flipud(im((R-padBR+1):R,:));
paddedImg((offset+1):(offset+R),(padC-padBR+1):padC) = fliplr(im(:,(C-padBR+1):C));
% now the corners
paddedImg(1:padTL,1:padTL) = fliplr(flipud(im(1:padTL,1:padTL)));
paddedImg(1:padTL,(padC-padBR+1):padC) = fliplr(flipud(im(1:padTL,(C-padBR+1):C)));
paddedImg((padR-padBR+1):padR,(padC-padBR+1):padC) = fliplr(flipud(im((R-padBR+1):R,(C-padBR+1):C)));
paddedImg((padR-padBR+1):padR,1:padTL) = fliplr(flipud(im((R-padBR+1):R,1:padTL)));

%paddedImg	% debug

% Now apply the spatial filtering
imOut = zeros(padR,padC);
for i = (offset+1):(offset+R)
   for j = (offset+1):(offset+C)
      imOut(i,j) = sum(sum(paddedImg((i-padTL):(i+padBR),(j-padTL):(j+padBR)) .* A));
   end
end

imOut = imOut((offset+1):(offset+R),(offset+1):(offset+C));

Sobel.m


function [gm,em] = sobel(im, T)
% [GM,EM] = SOBEL(IMG,T)
% Apply the Sobel gradient operator.  
% GM is the gradient magnitude
% EM is the edge map based on threshold T (T in range (0,1))
% The actual threshold value is determined such that 100*T% of pixels in gradient
% are kept as edges (T defaults to .10).

if nargin<2
   T = .10;
end

H1 =	[-1 0 1
   	-2 0 2
      -1 0 1];
   
H2 =	H1';

[R,C] = size(im);

N = 3;

% pad near the edges with the symmetric extension
if mod(N,2)
   padTL = (N-1)/2;	% left and top
   padBR = (N-1)/2;	% bottom and right
else
   padTL = N/2 - 1;
   padBR = N/2;
end
offset = padTL;
paddedIm = zeros(R+N-1, C+N-1);
[padR,padC] = size(paddedIm);
% copy the main image
paddedImg((offset+1):(offset+R), (offset+1):(offset+C)) = im;
% now copy up into the symmetric extension - top, left, bottom, right
% what a pain - is there a simpler way? (yes-sympad. thanks Lexing. but since
%   this works already, it's easier to leave it... :)
paddedImg(1:padTL,(offset+1):(offset+C)) = flipud(im(1:padTL,:));
paddedImg((offset+1):(offset+R),1:padTL) = fliplr(im(:,1:padTL));
paddedImg((padR-padBR+1):padR,(offset+1):(offset+C)) = flipud(im((R-padBR+1):R,:));
paddedImg((offset+1):(offset+R),(padC-padBR+1):padC) = fliplr(im(:,(C-padBR+1):C));
% now the corners
paddedImg(1:padTL,1:padTL) = fliplr(flipud(im(1:padTL,1:padTL)));
paddedImg(1:padTL,(padC-padBR+1):padC) = fliplr(flipud(im(1:padTL,(C-padBR+1):C)));
paddedImg((padR-padBR+1):padR,(padC-padBR+1):padC) = fliplr(flipud(im((R-padBR+1):R,(C-padBR+1):C)));
paddedImg((padR-padBR+1):padR,1:padTL) = fliplr(flipud(im((R-padBR+1):R,1:padTL)));

%paddedImg	% debug

% Apply horizontal gradient 
g1 = zeros(padR,padC);
for i = (offset+1):(offset+R)
   for j = (offset+1):(offset+C)
      g1(i,j) = paddedImg(i,j) - sum(sum(H1 .* paddedImg((i-padTL):(i+padBR),(j-padTL):(j+padBR))));
   end
end

% Apply vertical gradient 
g2 = zeros(padR,padC);
for i = (offset+1):(offset+R)
   for j = (offset+1):(offset+C)
      g1(i,j) = sum(sum(H2 .* paddedImg((i-padTL):(i+padBR),(j-padTL):(j+padBR))));
   end
end

% gradient magnitude:
gm = sqrt(g1.^2 + g2.^2);
% discard padded extensions
gm = gm((offset+1):(offset+R),(offset+1):(offset+C));

% normalize back to [0,1]
%gm = gm ./ max(max(gm));

% threshold to get the edge map
THRESH_PERCENT = 1-T;	% so 100*T% of pixels will be kept
NBINS = 1000;
% define threshold dynamically so that there are THRESH_PERCENT pixels in the edge map
[h,bins] = hist(gm(:),NBINS);
ch = cumsum(h);
% threshold cum pdf 
ch = find(ch > THRESH_PERCENT * sum(h));
t = bins(ch(1));
em = (gm > t);

edgemap.m


function em = edgemap(g,T)
% EM= EDGEMAP(G,T)
%   threshold the gradient G at 100*T% to get the edge map EM
% Adam Berenzweig, April 2001

% threshold to get the edge map
THRESH_PERCENT = 1-T;	% so 100*T% of pixels will be kept
NBINS = 1000;
% define threshold dynamically so that there are THRESH_PERCENT of pixels in the edge map
[h,bins] = hist(g(:),NBINS);	% abs cuz laplacian gives neg and pos values
ch = cumsum(h);
% threshold cum pdf 
ch = find(ch > THRESH_PERCENT * sum(h));
t = bins(ch(1));
em = (g > t);

laplace.m


function [g,em] = laplace(im,T,op)
% [L,EM] = LAPLACE(IMG,T,OP)
% Apply the Laplacian gradient operator.  
% L is the Laplacian gradient
% EM is the edge map based on threshold T (T in range (0,1))
% The actual threshold value is determined such that 100*T% of pixels in gradient
% are kept as edges (T defaults to .10).
% OP is which discrete laplacian to use:  (default is L2):
%  L1 = [ 0 -1 0 
%        -1 4 -1
%         0 -1 0 ];
%   
%  L2 = [ -1 -1 -1
%         -1 8 -1
%         -1 -1 -1 ];
%   
%  L3 = [ 1 -2 1
%        -2 4 -2
%         1 -2 1 ];

if nargin<2
   T = .10;
end

L1 = [ 0 -1 0 
   	-1 4 -1
      0 -1 0 ];
   
L2 = [ -1 -1 -1
		-1 8 -1
      -1 -1 -1 ];
   
L3 = [ 1 -2 1
		-2 4 -2
      1 -2 1 ];
   
% choose which laplacian from table 4:
if nargin < 3
   L = L2;
else
   Ls = {L1 L2 L3};
   L = Ls{op};
end
L

[R,C] = size(im);

N = 3;

% pad near the edges with the symmetric extension
if mod(N,2)
   padTL = (N-1)/2;	% left and top
   padBR = (N-1)/2;	% bottom and right
else
   padTL = N/2 - 1;
   padBR = N/2;
end
offset = padTL;
paddedIm = zeros(R+N-1, C+N-1);
[padR,padC] = size(paddedIm);
% copy the main image
paddedImg((offset+1):(offset+R), (offset+1):(offset+C)) = im;
% now copy up into the symmetric extension - top, left, bottom, right
% what a pain - is there a simpler way? (yes-sympad. thanks Lexing. but since
%   this works already, it's easier to leave it... :)
paddedImg(1:padTL,(offset+1):(offset+C)) = flipud(im(1:padTL,:));
paddedImg((offset+1):(offset+R),1:padTL) = fliplr(im(:,1:padTL));
paddedImg((padR-padBR+1):padR,(offset+1):(offset+C)) = flipud(im((R-padBR+1):R,:));
paddedImg((offset+1):(offset+R),(padC-padBR+1):padC) = fliplr(im(:,(C-padBR+1):C));
% now the corners
paddedImg(1:padTL,1:padTL) = fliplr(flipud(im(1:padTL,1:padTL)));
paddedImg(1:padTL,(padC-padBR+1):padC) = fliplr(flipud(im(1:padTL,(C-padBR+1):C)));
paddedImg((padR-padBR+1):padR,(padC-padBR+1):padC) = fliplr(flipud(im((R-padBR+1):R,(C-padBR+1):C)));
paddedImg((padR-padBR+1):padR,1:padTL) = fliplr(flipud(im((R-padBR+1):R,1:padTL)));

%paddedImg	% debug

% Apply laplacian
g = zeros(padR,padC);
for i = (offset+1):(offset+R)
   for j = (offset+1):(offset+C)
      g(i,j) = sum(sum(L .* paddedImg((i-padTL):(i+padBR),(j-padTL):(j+padBR))));
   end
end

% discard padded extensions
g = g((offset+1):(offset+R),(offset+1):(offset+C));

% normalize back to [0,1]
%g = g ./ max(max(g));

% threshold to get the edge map
THRESH_PERCENT = 1-T;	% so 100*T% of pixels will be kept
NBINS = 1000;
% de-mean to center the histogram
g = g - mean(g(:));
% define threshold dynamically so that there are THRESH_PERCENT pixels in the edge map
[h,bins] = hist(abs(g(:)),NBINS);	% abs cuz laplacian gives neg and pos values
ch = cumsum(h);
% threshold cum pdf 
ch = find(ch > THRESH_PERCENT * sum(h));
t = bins(ch(1));
em = (abs(g) > t);

dbledgemap.m


function em = dbledgemap(g,T)
% EM= DBLEDGEMAP(G,T)
%   threshold the (laplacian) gradient G at 100*T% to get the
%   (double) edge map EM
% Adam Berenzweig, April 2001

% threshold to get the edge map
THRESH_PERCENT = 1-T;	% so 100*T% of pixels will be kept
NBINS = 1000;

% de-mean to center the histogram
%g = g - mean(g(:));

% define threshold dynamically so that there are THRESH_PERCENT of pixels in the edge map
[h,bins] = hist(abs(g(:)),NBINS);	% abs cuz laplacian gives neg and pos values
ch = cumsum(h);
% threshold cum pdf 
ch = find(ch > THRESH_PERCENT * sum(h));
t = bins(ch(1));
em = (abs(g) > t);

Back