% These are commands you can use in order to test the different region
% based segmentation algorithms. Paste them in the Matlab command window
% and things should work.
% Test the merging algorithm
% Clear everything
clear all
close all
% Go to your image catalog and read your image (exchange with your own
% catalog)
i=double(imread('cameraman.tif'));
%i=double(rgb2gray(imread('mm3.tif')));
% Take a look at the image
figure
imshow(i,[min(min(i)) max(max(i))])
% Make initial suggestion for segmented image
tmp_seg=reshape(1:prod(size(i)),size(i));
% Test elementary region growing on this image. Try changing the value
% for the criterion and the connectivity. And obviously, try different
% images.
[seg,seg_arr,reg]=region_grow1(i,tmp_seg,25,8);
% Take a look at the result
map=rand(max(max(seg)),3);
figure
imshow(seg,[])
colormap(map)
function [seg,seg_arr,reg_info]=region_grow1(im,seg,criterion,connectivity)
% This function performs basic region growing on an image.
%
% The input arguments are the following:
% im: The image to be segmented
% seg: An initial labelling of the image, typically you would let
% this image be an image where every pixel has its own label
% criterion: The value of the similarity criterion that is used to decide
% between fusion of segments or not.
% connectivity: The connectivity of regions considered for fusion.
%
% The output arguments are the following:
%
% seg: The final segmented image
% seg_arr: A stack of all the segmented images from every
% iteration
% reg_info: A structure containing all the information about
% the different segments.
%
% Copyright Lars Aurdal, The Norwegian Computing Centre
% Get size of input image
rows=size(im,1);
cols=size(im,2);
% Here comes a slightly tricky part. In order to keep the computation
% time down we need to use s data structure to keep track of information
% about the regions. A MATLAB struct is what we need. Each element in
% the reg_info array of structures holds the size in pixels of the
% element, the mean of the pixels assigned to the element as well as the
% row and column coordinates of all the pixels assigned to the region,
% this is nice to have in order to be able to quickly adress all the
% pixels belonging to a single element.
% Start by initialising
reg_info=struct('size',num2cell(zeros(1,max(max(seg)))),... % Initially we don't know their size
'crt',num2cell(zeros(1,max(max(seg)))),... % ...nor their criterion value
'row_coord',[],... % Coordinates go here
'col_coord',[]); % and here
% Now loop over all the pixels in the initial segmented image and gather
% information about the existing regions
for m=1:rows
for n=1:cols
c_lab=seg(m,n);
reg_info(c_lab).size = reg_info(c_lab).size+1;
% The following calculation of criterion corresponds to a calculation
% of the segment means
reg_info(c_lab).crt = ((reg_info(c_lab).size-1)*reg_info(c_lab).crt+...
im(m,n))/reg_info(c_lab).size;
reg_info(c_lab).row_coord = [reg_info(c_lab).row_coord m];
reg_info(c_lab).col_coord = [reg_info(c_lab).col_coord n];
end
end
% Initially we assume that something will change in the first iteration
changed=1;
% Keep track of iterations in order to display this
iter=0;
% Now loop over the segmented image as long as something changes.
while(changed)
% Give user some feedback
disp(['Currently in iteration: ' num2str(iter+1)]);
% Check if something changes in this iteration
changed=0;
% Now loop
for m=1:rows
for n=1:cols
% Get the label of the current pixel considered
c_lab=seg(m,n);
% Get the coordinates of the current pixels neighbours. This
% can be done in 4 or 8 connectivity
if (connectivity==4)
ngb=get4ngb(rows,cols,m,n);
end
if (connectivity==8)
ngb=get8ngb(rows,cols,m,n);
end
% Loop over all these neighbours and check if they should
% receive the same label as the pixel we are currently
% examining
for k=1:size(ngb,2)
n_lab=seg(ngb(1,k),ngb(2,k));
if(n_lab~=c_lab)
if(abs(reg_info(n_lab).crt-reg_info(c_lab).crt)<criterion)
% First update the mean of the segment having the
% c_lab label
reg_info(c_lab).crt=(reg_info(c_lab).size*reg_info(c_lab).crt+...
reg_info(n_lab).size*reg_info(n_lab).crt)/...
(reg_info(c_lab).size+reg_info(n_lab).size);
% Then we are ready to update the size
reg_info(c_lab).size=reg_info(c_lab).size+reg_info(n_lab).size;
% Then give all pixels with n_lab the c_lab label
tmp_row_coord=reg_info(n_lab).row_coord;
tmp_col_coord=reg_info(n_lab).col_coord;
for p=1:reg_info(n_lab).size
seg(tmp_row_coord(p),tmp_col_coord(p))=c_lab;
end
% Finally update the coordinates of the pixels
% belonging to this segment
reg_info(c_lab).row_coord=[reg_info(c_lab).row_coord reg_info(n_lab).row_coord];
reg_info(c_lab).col_coord=[reg_info(c_lab).col_coord reg_info(n_lab).col_coord];
% ...and mark it for deletion
reg_info(n_lab).size=0;
% If we have reached this point then we know that
% something changed in this iteration
changed=1;
end
end
end
end
end
iter=iter+1;
% We will store the temporary results in seg_arr for demonstration
% purposes, the first layer equals the initial segmentation
seg_arr(:,:,iter)=seg;
end
% Clean up the reg_info structure
index=1;
for m=1:max(max(seg))
if (reg_info(m).size>0)
tmp_reg_info(index)=reg_info(m);
index=index+1;
end
end
reg_info=tmp_reg_info;
function ngb = get4ngb(rows,cols,x,y)
% function ngb = get4ngb(rows,cols,x,y)
% This function returns the cooridinates of the 8-neighbours
% of an element in a matrix. The values are returned in the
% list ngb. If the neighbour does not exist,- that is (x,y)
% corredsponds to an edge or corner of the array,
% the list of neighbours contains only those coordinates corresponding
% to real neighbours.
% Copyright Lars Aurdal/FFIE.
% Handle left edge.
% Handle upper left corner.
if ((x == 1) & (y == 1))
ngb(1:2,1:1) = [2 1]';
ngb(1:2,2:2) = [1 2]';
% Handle lower left corner.
elseif ((x == rows) & (y == 1))
ngb(1:2,1:1) = [rows 2]';
ngb(1:2,2:2) = [(rows - 1) 1]';
% Handle left edge in general.
elseif (y == 1)
ngb(1:2,1:1) = [(x+1) 1]';
ngb(1:2,2:2) = [x 2]';
ngb(1:2,3:3) = [(x-1) 1]';
% Handle right edge.
% Handle upper right corner.
elseif ((x == 1) & (y == cols))
ngb(1:2,1:1) = [1 (cols-1)]';
ngb(1:2,2:2) = [2 cols]';
% Handle lower right corner.
elseif ((x == rows) & (y == cols))
ngb(1:2,1:1) = [rows (cols-1)]';
ngb(1:2,2:2) = [(rows-1) cols]';
% Handle right edge in general.
elseif (y == cols)
ngb(1:2,1:1) = [(x+1) cols]';
ngb(1:2,2:2) = [x (cols-1)]';
ngb(1:2,3:3) = [(x-1) cols]';
% Handle top line.
elseif (x == 1)
ngb(1:2,1:1) = [1 (y-1)]';
ngb(1:2,2:2) = [2 y]';
ngb(1:2,3:3) = [1 (y+1)]';
% Handle bottom line.
elseif (x == rows)
ngb(1:2,1:1) = [rows (y-1)]';
ngb(1:2,2:2) = [(rows-1) y]';
ngb(1:2,3:3) = [rows (y+1)]';
% Handle general case
else
ngb(1:2,1:1) = [(x-1) y]';
ngb(1:2,2:2) = [x (y-1)]';
ngb(1:2,3:3) = [(x+1) y]';
ngb(1:2,4:4) = [x (y+1)]';
end
function ngb = get8ngb(rows,cols,x,y)
% function ngb = get8ngb(rows,cols,x,y)
% This function returns the cooridinates of the 8-neighbours of an element
% in a matrix. The values are returned in the list neighbours. If the neighbour
% does not exist,- that is center corredsponds to an edge or corner of the array,
% the list of neighbours contains only those coordinates corresponding to real
% neighbours.
% Copyright Lars Aurdal/FFIE.
% Handle left edge.
% Handle upper left corner.
if ((x == 1) & (y == 1))
ngb(1:2,1:1) = [2 1]';
ngb(1:2,2:2) = [2 2]';
ngb(1:2,3:3) = [1 2]';
% Handle lower left corner.
elseif ((x == rows) & (y == 1))
ngb(1:2,1:1) = [rows 2]';
ngb(1:2,2:2) = [(rows - 1) 2]';
ngb(1:2,3:3) = [(rows - 1) 1]';
% Handle left edge in general.
elseif (y == 1)
ngb(1:2,1:1) = [(x+1) 1]';
ngb(1:2,2:2) = [(x+1) 2]';
ngb(1:2,3:3) = [x 2]';
ngb(1:2,4:4) = [(x-1) 2]';
ngb(1:2,5:5) = [(x-1) 1]';
% Handle right edge.
% Handle upper right corner.
elseif ((x == 1) & (y == cols))
ngb(1:2,1:1) = [1 (cols-1)]';
ngb(1:2,2:2) = [2 (cols-1)]';
ngb(1:2,3:3) = [2 cols]';
% Handle lower right corner.
elseif ((x == rows) & (y == cols))
ngb(1:2,1:1) = [rows (cols-1)]';
ngb(1:2,2:2) = [(rows-1) (cols-1)]';
ngb(1:2,3:3) = [(rows-1) cols]';
% Handle right edge in general.
elseif (y == cols)
ngb(1:2,1:1) = [(x+1) cols]';
ngb(1:2,2:2) = [(x+1) (cols-1)]';
ngb(1:2,3:3) = [x (cols-1)]';
ngb(1:2,4:4) = [(x-1) (cols-1)]';
ngb(1:2,5:5) = [(x-1) cols]';
% Handle top line.
elseif (x == 1)
ngb(1:2,1:1) = [1 (y-1)]';
ngb(1:2,2:2) = [2 (y-1)]';
ngb(1:2,3:3) = [2 y]';
ngb(1:2,4:4) = [2 (y+1)]';
ngb(1:2,5:5) = [1 (y+1)]';
% Handle bottom line.
elseif (x == rows)
ngb(1:2,1:1) = [rows (y-1)]';
ngb(1:2,2:2) = [(rows-1) (y-1)]';
ngb(1:2,3:3) = [(rows-1) y]';
ngb(1:2,4:4) = [(rows-1) (y+1)]';
ngb(1:2,5:5) = [rows (y+1)]';
% Handle general case
else
ngb(1:2,1:1) = [(x-1) y]';
ngb(1:2,2:2) = [(x-1) (y-1)]';
ngb(1:2,3:3) = [x (y-1)]';
ngb(1:2,4:4) = [(x+1) (y-1)]';
ngb(1:2,5:5) = [(x+1) y]';
ngb(1:2,6:6) = [(x+1) (y+1)]';
ngb(1:2,7:7) = [x (y+1)]';
ngb(1:2,8:8) = [(x-1) (y+1)]';
end
代码
