|
Vision Project - CS223B - Winter 97 Eric Frew, Andreas Huster, Edward LeMaster Matlab Code
Convolution Kernel Generationfunction kernel=sdg(sigma) % % kernel=sdg(sigma) % % Determines a derivative of gaussian convolution kernel from sigma % 'sigma' must be an integer % Width of the kernel is 6*sigma-1 % % Author: EAL 2-97 a=[-3*sigma+1:1:3*sigma-1]; kernel=-a.*exp(-.5*(a/sigma).^2)/(sigma^2); Histogram Local Minima
function index=local_min(histogram)
%
% index=locl_min(histogram)
%
% Determines the local minima of a sdg(histogram) vector
% Convolution with sdg kernel must take place previously
% Uses zero-crossings of the vector
% index is a vector giving the indices of all minima
%
% Author: EAL 2-97
histogram=sign(histogram);
prev=1;
for i=1:length(histogram),
current=histogram(i);
if current==1,
if prev 'less than or =' 0,
index=[index;i];
end
end
prev=current;
end
function index2=min_thresh_loc(histogram,index,percent,pixels)
%
% index2=min_thresh_loc(histogram,index,percent)
%
% Determines which local minima are suitable for threshold locations.
% (Thresholds are not added unless 'bins' are or a certain size.)
% 'histogram' is the original (unfiltered) histogram
% 'index' is the vector of all local minima
% 'percent' is the fraction of pixels required between threshold locations
% 'pixels' is the total number of pixels in the image
%
% Author: EAL 2-97
min_bin=percent*pixels;
prev_index=0;
for i=1:length(index),
item=index(i);
s=sum(histogram((prev_index+1):item));
str=sprintf('sum(%d)=%d',item,s);
disp(str);
if s >= min_bin,
index2=[index2; item];
prev_index=item;
end
end
Histogram Algorithm #1
% Loads in pgm images and determines automatic thresholding.
% Outputs histograms and thresholded images
%
% Author: EAL 2-97
%%% Get Image Number %%%
image_number=input('Please input the image number >');
%%% Load Pictures %%%
str=sprintf('l%d.pgm',image_number);
l_image=pgmread(str);
str=sprintf('u%d.pgm',image_number);
u_image=pgmread(str);
str=sprintf('v%d.pgm',image_number);
v_image=pgmread(str);
[m,n]=size(l_image);
[m2,n2]=size(u_image);
%%% Pseudo-Gaussian Smoothing %%%
kernel=[1/16 1/8 1/16; 1/8 1/4 1/8; 1/16 1/8 1/16];
kernel2=[0 1/8 0;1/8 1/2 1/8; 0 1/8 0];
l_image2=conv2(l_image,kernel,'same');
u_image2=conv2(u_image,kernel,'same');
v_image2=conv2(v_image,kernel2,'same');
%%% Convert Matrices to Vectors %%%
l_vec=reshape(l_image,[1 m*n]); % Used for debugging
u_vec=reshape(u_image,[1 m2*n2]); % Used for debugging
v_vec=reshape(v_image,[1 m2*n2]); % Used for debugging
l_vec2=reshape(l_image2,[1 m*n]);
u_vec2=reshape(u_image2,[1 m2*n2]);
v_vec2=reshape(v_image2,[1 m2*n2]);
%%% Generate Histograms %%%
bin=[0:1:255];
[nl,xl]=hist(l_vec,bin); % Used for debugging
[nu,xu]=hist(u_vec,bin); % Used for debugging
[nv,xv]=hist(v_vec,bin); % Used for debugging
[nl2,xl2]=hist(l_vec2,bin);
[nu2,xu2]=hist(u_vec2,bin);
[nv2,xv2]=hist(v_vec2,bin);
%%% Smooth Histograms %%%
h_kernel=sdg(5);
nl2_smooth=conv2(nl2,h_kernel,'same');
nu2_smooth=conv2(nu2,h_kernel,'same');
nv2_smooth=conv2(nv2,h_kernel,'same');
%%% Determine U Thresholds %%%
indexl=local_min(nl2_smooth);
index2l=min_thresh_loc(nl2,indexl,0.20,length(l_vec2));
dummyl=zeros(size(index2l));
indexu=local_min(nu2_smooth);
index2u=min_thresh_loc(nu2,indexu,0.20,length(u_vec2));
dummyu=zeros(size(index2u));
indexv=local_min(nv2_smooth);
index2v=min_thresh_loc(nv2,indexv,0.20,length(v_vec2));
dummyv=zeros(size(index2v));
%%% Plot Histograms %%%
clg;
str=sprintf('Image %d Histograms - Algorithm #1',image_number);
subplot(311); hold on; title(str);
stairs(xl,nl); plot(xl2,nl2,'r'); plot(xl2,nl2_smooth,'b');
plot(index2l, dummyl,'x'); ylabel('L-image');
subplot(312); hold on;
stairs(xu,nu); plot(xu2,nu2,'r'); plot(xu2,nu2_smooth,'b');
plot(index2u, dummyu,'x'); ylabel('U-image');
subplot(313); hold on;
stairs(xv,nv); plot(xv2,nv2,'r'); plot(xv2,nv2_smooth,'b');
plot(index2v, dummyv,'x'); ylabel('V-image');
%%% Auto-Threshold and Save %%%
if length(index2l)>0,
k=length(index2l);
temp=im2bw(l_image2,index2l(k));
str=sprintf('l%dt%d.pgm',image_number,k);
pgmwrite(temp,str);
end
if length(index2u)>0,
k=length(index2u);
temp=im2bw(u_image2,index2u(k));
str=sprintf('u%dt%d.pgm',image_number,k);
pgmwrite(temp,str);
end
if length(index2v)>0,
k=length(index2v);
temp=im2bw(v_image2,index2v(k));
str=sprintf('v%dt%d.pgm',image_number,k);
pgmwrite(temp,str);
end
Histogram Algorithm #2Code to load images and call algorithm:
%CS223B, Project, Pixel Selection based on a histogram of pixels with a high
% gradient
clear;
%Get image number
image_number=input('Please input the image number >');
%path for images
datapath = sprintf('../data');
resultpath = sprintf('../results');
%load pictures
str=sprintf('%s/u%d.pgm', datapath, image_number);
u_image=pgmread(str);
str=sprintf('%s/v%d.pgm', datapath, image_number);
v_image=pgmread(str);
%detect objects
[u_result, u_grad_th] = process_image(u_image, 1);
[v_result, v_grad_th] = process_image(v_image, 2);
%save result images
str=sprintf('%s/u%d.ah.pgm', resultpath, image_number);
pgmwrite(u_result, str);
str=sprintf('%s/v%d.ah.pgm', resultpath, image_number);
pgmwrite(v_result, str);
%construct and save the composite image
[a b]=size(u_result);
str=sprintf('%s/c%d.ah.pgm', resultpath, image_number);
c = [u_result ones(a, 1) v_result;
ones(1,2*b+1);
u_grad_th ones(a, 1) v_grad_th];
pgmwrite(c, str);
Actual Algorithm:
function [result,grad_th]=process_image(image, object_type)
%PROCESS_IMAGE This function attempts to extract objects from
% the image using the following procedure:
% 1. smooth the image
% 2. calculate image gradient
% 3. calculate histogram of gradient
% 4. calculate gradient cutoff
% 5. determine set of pixels with high gradients
% 6. calculate image threshold mean of this set
% 7. Determine validity of threshold
% 8. Apply threshold to image
%
% Object Types are:
% 1. Barrel
% 2. Disk
% Andreas Huster Feb 26, 1997
%%% MAGIC NUMBERS
%percentile cutoff for gradients
if (object_type == 1)
grad_percentile = 0.03;
obj_str = sprintf('barrel');
else
grad_percentile = 0.01;
obj_str = sprintf('disk');
end;
ker_size = 2; %kernel size for smoothing
NGradBins = 30; %number of bins in gradiant histogram
%image size
[m n]=size(image);
%%%% SMOOTHING
gauss_kernel = [0 1/32 1/16 1/32 0;
1/32 1/32 1/16 1/32 1/32;
1/16 1/16 1/8 1/16 1/16;
1/32 1/32 1/16 1/32 1/32;
0 1/32 1/16 1/32 0];
%smooth
smooth = conv2(image, gauss_kernel, 'same');
%%%% GRADIENT
grad_kernel = [ 0 1/2*i 0;
1/2 0 -1/2;
0 -1/2*i 0];
grad_uncropped = abs(conv2(smooth, grad_kernel, 'same'));
%Crop the edge of the image where the gradient is poorly defined
crop_size = ker_size*3;
grad = grad_uncropped((1+crop_size):(m-crop_size), ...
(1+crop_size):(n-crop_size));
%crop the smoothed image as well so that index values work out
csmooth = smooth((1+crop_size):(m-crop_size), ...
(1+crop_size):(n-crop_size));
mc = m-2*crop_size;
nc = n-2*crop_size;
%%%% GRADIENT HISTOGRAM
%generate a vector of gradient values for ease of use in histograms
grad_vec = reshape(grad, [mc*nc 1]);
%histogram the gradient values
[N,X]=hist(grad_vec,NGradBins);
%find the gradient cutoff for the nth percentile
Nmin = grad_percentile*length(grad_vec);
cutoff_bin = NGradBins;
Nover = N(NGradBins);
%loop until we contain enough bins
while (Nover < Nmin)
cutoff_bin = cutoff_bin - 1;
Nover = Nover + N(cutoff_bin);
end;
grad_cutoff = X(cutoff_bin);
%%%% THRESHOLD CALCULATION
%find the high-gradient pixels in cropped image
[I,J]=find(grad > grad_cutoff);
Index = mc*(J-1)+I;
%find corresponding pixels in the smoothed image
high_grad_pixels = csmooth(Index);
%compute mean of edge pixels
image_cutoff = mean(high_grad_pixels);
%%%% DETECTION LOGIC
mean_image = mean(mean(image));
if (image_cutoff < mean_image*1.1)
image_cutoff = 300;
str=sprintf('No %ss detected.', obj_str);
else
str=sprintf('%s: cutoff = %4.0f',obj_str,image_cutoff);
end;
disp(str);
%%%% THRESHOLDING
result = im2bw(smooth,image_cutoff);
grad_th = im2bw(grad_uncropped, grad_cutoff);
Percentile Threshold Determinationfunction thresh=ptile_thresh_find(hist_n,hist_x,total_pixels,ptile) % This function returns a threshold such that ptile percent of the total % number of pixels falls above the threshold. % Author: Eric Frew 2-97 num_pixels=total_pixels*ptile; i=length(hist_x); above_thresh=0; while(above_thresh<=num_pixels) add_pixels=hist_n(i); above_thresh=above_thresh+add_pixels; i=i-1; end thresh=hist_x(i); Percentile Algorithm #3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eric Frew
% 2-26-96
% loads in PGM versions of LUV images and then computes histograms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
flag=0;
%%% parameters %%%
u_ptile=.08;
v_ptile=.0025;
image_number=1;
im_nums=[10 11 14 15 16 17 18 19 21 22 24 25 3 5];
for num=1:length(im_nums)%%%
%%% Get image number %%%
image_number=im_nums(num)
%image_number=input('Please input the image number >');
%%% load pictures %%%
str=sprintf('l%d.pgm',image_number);
l_image=pgmread(str);
str=sprintf('u%d.pgm',image_number);
u_image=pgmread(str);
str=sprintf('v%d.pgm',image_number);
v_image=pgmread(str);
[m,n]=size(l_image);
[m2,n2]=size(u_image);
%%% Gaussian Smoothing %%%
kernel=[1/16 1/8 1/16; 1/8 1/4 1/8; 1/16 1/8 1/16];
kernel2=[0 1/8 0;1/8 1/2 1/8; 0 1/8 0];
l_image2=conv2(l_image,kernel,'same');
u_image2=conv2(u_image,kernel,'same');
v_image2=conv2(v_image,kernel2,'same');
%%% Convert Matrices to Vectors %%%
l_vec2=reshape(l_image2,[1 m*n]);
u_vec2=reshape(u_image2,[1 m2*n2]);
v_vec2=reshape(v_image2,[1 m2*n2]);
%%% Generate Histograms %%%
bin=[0:1:255];
[nl2,xl2]=hist(l_vec2,bin);
[nu2,xu2]=hist(u_vec2,bin);
[nv2,xv2]=hist(v_vec2,bin);
%%% Determine Threshold %%%
thresh_u2=ptile_thresh_find(nu2,xu2,m2*n2,u_ptile);
thresh_v2=ptile_thresh_find(nv2,xv2,m2*n2,v_ptile);
if flag==1
%%% Threshold images %%%
u_im_thresh=im2bw(u_image2,thresh_u2);
v_im_thresh=im2bw(v_image2,thresh_v2);
%%% Smooth Threshold (eliminate non-barrels) %%%
u_im_smooth=conv2(255*u_im_thresh,kernel,'same');
v_im_smooth=conv2(255*v_im_thresh,kernel2,'same');
%%% Smooth Histograms %%%
h_kernel=sdg(5);
nl2_smooth=conv2(nl2,h_kernel,'same');
nu2_smooth=conv2(nu2,h_kernel,'same');
nv2_smooth=conv2(nv2,h_kernel,'same');
%%% Plot Histograms %%%
figure(1)
clg;
str=sprintf('Image %d Histograms - Algorithm #1',image_number);
subplot(311); hold on; title(str);
plot(xl,nl,'g'); plot(xl2,nl2,'r');
ylabel('L-image');
subplot(312); hold on;
plot(xu,nu,'g'); plot(xu2,nu2,'r');
plot(thresh_u2, nu2(thresh_u2+1),'mx');
ylabel('U-image');
subplot(313); hold on;
plot(xv,nv,'g'); plot(xv2,nv2,'r');
plot(thresh_v2, nv2(thresh_v2-1),'mx');
ylabel('V-image');
legend('Original Image','Smoothed Image');
end %flag
%%% Auto-Threshold and Save %%%
if length(thresh_u2)>1,
for k=1:length(index2l),
temp=im2bw(l_image2,index2l(k));
str=sprintf('l%dt%d.pgm',image_number,k);
pgmwrite(temp,str);
end
end
if length(thresh_u2)>0,
temp=im2bw(u_image2,thresh_u2);
str=sprintf('u%dt%d.pgm',image_number,1);
pgmwrite(temp,str);
end
if length(thresh_v2)>0,
temp=im2bw(v_image2,thresh_v2);
str=sprintf('v%dt%d.pgm',image_number,1);
pgmwrite(temp,str);
end
end % for num=...
SegmentationCode to load images and call algorithm:
%CS223B, Project
%
%Segments the regions in the barrel/disk b/w images and figures
%the centroid marking it with a + or x.
clear;
%Get image number
image_number=input('Please input the image number: ');
%path for images
resultpath = sprintf('../results');
%load pictures
str=sprintf('%s/u%d.ah.pgm', resultpath, image_number);
u_bw=pgmread(str);
str=sprintf('%s/v%d.ah.pgm', resultpath, image_number);
v_bw=pgmread(str);
segs = zeros(size(u_bw));
%detect objects
disp(' ');
disp('Object row column size');
disp('----------------------------------');
segs = find_objects(u_bw, segs, 1);
segs = find_objects(v_bw, segs, 2);
%save result
str=sprintf('%s/seg%d.pgm', resultpath, image_number);
pgmwrite(segs, str);
Actual Algorithm:
function [segs]=find_objects(bw, prev_segs, object_type)
%PROCESS_IMAGE This function segments the b/w image using
% the following procedure:
% 1. while there are white pixels that haven't been
% accounted for
% 2. initialize open list with one of the
% unaccounted pixels
% 3. while there are pixels in the open list
% 4. take the next pixel in the list and
% enumerate its neighbors
% 5. add neighbors which are white and not
% already in the open list to the list
% 6. make the white neighbors black
% 7. all pixels that were in the open list make
% up one region
% 8. calculate centroid of region
%
% Object Types are:
% 1. Barrel
% 2. Disk
% Andreas Huster Mar 9, 1997
%Marker for objects
if (object_type == 1)
% an x
marker = [1 1 0 0 0 0 0 1 1;
1 1 1 0 0 0 1 1 1;
0 1 1 1 0 1 1 1 0;
0 0 1 1 1 1 1 0 0;
0 0 0 1 1 1 0 0 0;
0 0 1 1 1 1 1 0 0;
0 1 1 1 0 1 1 1 0;
1 1 1 0 0 0 1 1 1;
1 1 0 0 0 0 0 1 1];
mSize = 4;
else
% a plus
marker = [0 0 0 0 1 1 1 0 0 0 0;
0 0 0 0 1 1 1 0 0 0 0;
0 0 0 0 1 1 1 0 0 0 0;
0 0 0 0 1 1 1 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 1 1 1 0 0 0 0;
0 0 0 0 1 1 1 0 0 0 0;
0 0 0 0 1 1 1 0 0 0 0;
0 0 0 0 1 1 1 0 0 0 0];
mSize = 5;
end;
%add to the previous markers
segs = prev_segs;
%these vectors define the relative coordinates of the neighborhood
%of a pixel
NI = [0, 1, 0, -1];
NJ = [1, 0, -1, 0];
numN = 4;
%find the total size of the image
[iMax,jMax] = size(bw);
%find a pixel that is part of one of the regions
[I,J] = find(bw > 0);
%loop for as long as there are still new regions
while (length(I) > 0)
%initialize open list with one of the pixels
openI = I(1);
openJ = J(1);
cur = 1;
%loop until the entire open list is processed
while (cur <= length(openI))
%take the next point in the open list
i = openI(cur);
j = openJ(cur);
%enumerate all of its neighbors
for n = 1:numN
in = i + NI(n);
jn = j + NJ(n);
%make sure the neighbor pixel is in bounds
if ((in > 0) & (in <= iMax) & (jn > 0) & (jn <= jMax))
%check whether the pixel is part of the region
if (bw(in,jn) > 0)
%if yes, add it to the open list
openI = [openI; in];
openJ = [openJ; jn];
%and make it black in the image
bw(in,jn) = 0;
end;
end;
end;
%look at next pixel in the list
cur = cur + 1;
end;
%now we have a region
obj_size = length(openI);
center_i = mean(openI);
center_j = mean(openJ);
%display result
if (object_type == 1)
obj_str = sprintf('Barrel ');
else
obj_str = sprintf('Disk ');
end;
str = sprintf('%s%5.2f %5.2f %5d', ...
obj_str, center_i, center_j, obj_size);
disp(str);
%filter out barrels that are too small
if (~((object_type == 1) & (obj_size < 100)))
center_i = round(center_i);
center_j = round(center_j);
%place a marker on the segs image
r1 = (center_i - mSize): (center_i + mSize);
r2 = (center_j - mSize): (center_j + mSize);
segs(r1,r2) = segs(r1,r2) | marker;
end;
%look for any more unprocessed pixels
[I,J] = find(bw > 0);
end;
Return to the main page. Copyright © 1997 by Eric Frew, Andreas Huster, & Edward LeMaster. |