Homework 4: Matlab Code

• spatial average
• gradient operator
• Unsharp masking
• Threshold slicing
• Threshold slice with stretching
• Median Filtering
• Image normalization

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?
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));      

gradientfilt.m

function imOut = gradientfilt(im)
% B = GRADIENTFILT(A)		Gradient filter.
% The gradient uses the discrete laplacian;
% 	1/4 * [0 1 0
%			 1 0 1
%			 0 1 0]

% define discrete Laplacian for gradient operation. 
L = (1/4).*	[0 1 0;
   			 1 0 1;
             0 1 0];

[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?
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 gradient
imOut = zeros(padR,padC);
for i = (offset+1):(offset+R)
   for j = (offset+1):(offset+C)
      imOut(i,j) = paddedImg(i,j) - sum(sum(L .* paddedImg((i-padTL):(i+padBR),(j-padTL):(j+padBR))));
   end
end

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

function imOut = unsharp(im, G, lambda)
% O = UNSHARP(IM,G,L)		Unsharp masking
% Performs unsharp masking as follows:
% 		O = IM + L*G
% where G is the gradient performed by GRADIENTFILT(IM)
% O is returned normalized to (0,1)

%G = gradientfilt(im);
imOut = im + (lambda * G);

imOut = imnorm(imOut);

slice255.m

function imOut = slice255(im, a, b)
% O = SLICE(I, A, B)
% maps gray levels < A to 0, > B to 255, and preserves levels
% between A and B.  
% Using 255 gray levels, so 0255 | a>b)
   error('must have 0 <= A <= B <= 255');
end

mask = (im >=a & im <= b);		% 1's where we want slice.
imOut= im .* mask;
imOut = imOut + ~mask;

histstretch255.m

function imOut = histstretch255(im,a,b);
% OUT = HISTSTRETCH255(IMG,A,B)		Histogram stretching.  
% maps the grayscale range (A,B) of the input image to
% (0,255) in the output image.  values less than A are sent
% to 0, and values greater than B are sent to 255.

if (a < 0 | b>255 | a>b)
   error('must have 0 <= A <= B <= 255');
end

imOut = zeros(size(im));

% first mask range < a to 0
mask = (im < a);
imOut(mask) = 0;
% now range > b to 1
mask = (im > b);
imOut(mask) = 255;

% now expand range (a,b)
mask = (im >= a & im <= b);
d = (b-a);
imOut(mask) = floor((im(mask) - a) ./ d);

medianfilt.m

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

[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?
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 median filtering
imOut = zeros(padR,padC);
for i = (offset+1):(offset+R)
   for j = (offset+1):(offset+C)
      A = paddedImg((i-padTL):(i+padBR),(j-padTL):(j+padBR));
      imOut(i,j) = median(A(:));
   end
end

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

imnorm.m

function imOut = imnorm(im)
% O = IMNORM(I)		Image normalization
% imwrite doesnt seem to do any normalizing, so i do it here.
% normalizes to the range (0,1).

imOut = im;
% normalize
imMin = min(imOut(:));
if (imMin < 0)
   imOut = imOut + abs(imMin);
end
imOut = imOut ./ max(imOut(:));