Code covered by the BSD License  

Highlights from
Image segmentation using statistical region merging

image thumbnail

Image segmentation using statistical region merging

by

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);

% Compute image gradient
[Ix,Iy]=srm_imgGrad(image(:,:,:));
Ix=max(abs(Ix),[],3);
Iy=max(abs(Iy),[],3);
normgradient=sqrt(Ix.^2+Iy.^2);

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;
    bgradient = sparse(srm_boundarygradient(labels, nlabels, normgradient));
    bgradient = bgradient - tril(bgradient);
    idx=find(bgradient>0);
    [~,index]=sort(bgradient(idx));
    n_pairs=numel(idx);
    [xlabels,ylabels]=ind2sub([nlabels,nlabels],idx);
    pairs1=clusterlist(xlabels);
    pairs2=clusterlist(ylabels);
end










Contact us