区域增长

代码

% 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