区域增长

代码

% 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