Code covered by the BSD License

# Image segmentation using statistical region merging

### Sylvain Boltz (view profile)

20 Oct 2009 (Updated )

Very simple and fast image segmentation code using statistical region merging.

[maps,images]=srm(image,Qlevels)
```% Statistical Region Merging
%
% Nock, Richard and Nielsen, Frank 2004. Statistical Region Merging. IEEE Trans. Pattern Anal. Mach. Intell. 26, 11 (Nov. 2004), 1452-1458.
% DOI= http://dx.doi.org/10.1109/TPAMI.2004.110

%Segmentation parameter Q; Q small few segments, Q large may segments

function [maps,images]=srm(image,Qlevels)

% Smoothing the image, comment this line if you work on clean or synthetic images
h=fspecial('gaussian',[3 3],1);
image=imfilter(image,h,'symmetric');

smallest_region_allowed=10;

size_image=size(image);
n_pixels=size_image(1)*size_image(2);

Ix=max(abs(Ix),[],3);
Iy=max(abs(Iy),[],3);

Ix(:,end)=[];
Iy(end,:)=[];

[~,index]=sort(abs([Iy(:);Ix(:)]));

n_levels=numel(Qlevels);
maps=cell(n_levels,1);
images=cell(n_levels,1);
im_final=zeros(size_image);

map=reshape(1:n_pixels,size_image(1:2));
% gaps=zeros(size(map)); % For future release
treerank=zeros(size_image(1:2));

size_segments=ones(size_image(1:2));
image_seg=image;

%Building pairs
n_pairs=numel(index);
idx2=reshape(map(:,1:end-1),[],1);
idx1=reshape(map(1:end-1,:),[],1);

pairs1=[ idx1;idx2 ];
pairs2=[ idx1+1;idx2+size_image(1) ];

for Q=Qlevels
iter=find(Q==Qlevels);

for i=1:n_pairs
C1=pairs1(index(i));
C2=pairs2(index(i));

%Union-Find structure, here are the finds, average complexity O(1)
while (map(C1)~=C1 ); C1=map(C1); end
while (map(C2)~=C2 ); C2=map(C2); end

% Compute the predicate, region merging test
g=256;
logdelta=2*log(6*n_pixels);

dR=(image_seg(C1)-image_seg(C2))^2;
dG=(image_seg(C1+n_pixels)-image_seg(C2+n_pixels))^2;
dB=(image_seg(C1+2*n_pixels)-image_seg(C2+2*n_pixels))^2;

logreg1 = min(g,size_segments(C1))*log(1.0+size_segments(C1));
logreg2 = min(g,size_segments(C2))*log(1.0+size_segments(C2));

dev1=((g*g)/(2.0*Q*size_segments(C1)))*(logreg1 + logdelta);
dev2=((g*g)/(2.0*Q*size_segments(C2)))*(logreg2 + logdelta);

dev=dev1+dev2;

predicat=( (dR<dev) && (dG<dev) && (dB<dev) );

if (((C1~=C2)&&predicat) ||  xor(size_segments(C1)<=smallest_region_allowed, size_segments(C2)<=smallest_region_allowed))
% Find the new root for both regions
if treerank(C1) > treerank(C2)
map(C2) = C1; reg=C1;
elseif treerank(C1) < treerank(C2)
map(C1) = C2; reg=C2;
elseif C1 ~= C2
map(C2) = C1; reg=C1;
treerank(C1) = treerank(C1) + 1;
end

if C1~=C2
% Merge regions
nreg=size_segments(C1)+size_segments(C2);
image_seg(reg)=(size_segments(C1)*image_seg(C1)+size_segments(C2)*image_seg(C2))/nreg;
image_seg(reg+n_pixels)=(size_segments(C1)*image_seg(C1+n_pixels)+size_segments(C2)*image_seg(C2+n_pixels))/nreg;
image_seg(reg+2*n_pixels)=(size_segments(C1)*image_seg(C1+2*n_pixels)+size_segments(C2)*image_seg(C2+2*n_pixels))/nreg;
size_segments(reg)=nreg;
end
end
end

% Done, building two result figures, figure 1 is the segmentation map,
% figure 2 is the segmentation map with the average color in each segment

while 1
map_ = map(map) ;
if isequal(map_,map) ; break ; end
map = map_ ;
end

for i=1:3
im_final(:,:,i)=image_seg(map+(i-1)*n_pixels);
end
images{iter}=im_final;

[clusterlist,~,labels] = unique(map)  ;
labels=reshape(labels,size(map));
nlabels=numel(clusterlist);
maps{iter}=map;
n_pairs=numel(idx);
[xlabels,ylabels]=ind2sub([nlabels,nlabels],idx);
pairs1=clusterlist(xlabels);
pairs2=clusterlist(ylabels);
end

```