% SLIC Simple Linear Iterative Clustering SuperPixels
%
% Implementation of Achanta, Shaji, Smith, Lucchi, Fua and Susstrunk's
% SLIC Superpixels
%
% Usage: [l, Am, Sp, d] = slic(im, k, m, seRadius, colopt, mw)
%
% Arguments: im - Image to be segmented.
% k - Number of desired superpixels. Note that this is nominal
% the actual number of superpixels generated will generally
% be a bit larger, especially if parameter m is small.
% m - Weighting factor between colour and spatial
% differences. Values from about 5 to 40 are useful. Use a
% large value to enforce superpixels with more regular and
% smoother shapes. Try a value of 10 to start with.
% seRadius - Regions morphologically smaller than this are merged with
% adjacent regions. Try a value of 1 or 1.5. Use 0 to
% disable.
% colopt - String 'mean' or 'median' indicating how the cluster
% colour centre should be computed. Defaults to 'mean'
% mw - Optional median filtering window size. Image compression
% can result in noticeable artifacts in the a*b* components
% of the image. Median filtering can reduce this. mw can be
% a single value in which case the same median filtering is
% applied to each L* a* and b* components. Alternatively it
% can be a 2-vector where mw(1) specifies the median
% filtering window to be applied to L* and mw(2) is the
% median filtering window to be applied to a* and b*.
%
% Returns: l - Labeled image of superpixels. Labels range from 1 to k.
% Am - Adjacency matrix of segments. Am(i, j) indicates whether
% segments labeled i and j are connected/adjacent
% Sp - Superpixel attribute structure array with fields:
% L - Mean L* value
% a - Mean a* value
% b - Mean b* value
% r - Mean row value
% c - Mean column value
% stdL - Standard deviation of L*
% stda - Standard deviation of a*
% stdb - Standard deviation of b*
% N - Number of pixels
% edges - List of edge numbers that bound each
% superpixel. This field is allocated, but not set,
% by SLIC. Use SPEDGES for this.
% d - Distance image giving the distance each pixel is from its
% associated superpixel centre.
%
% It is suggested that use of this function is followed by SPDBSCAN to perform a
% DBSCAN clustering of superpixels. This results in a simple and fast
% segmentation of an image.
%
% Minor variations from the original algorithm as defined in Achanta et al's
% paper:
%
% - SuperPixel centres are initialised on a hexagonal grid rather than a square
% one. This results in a segmentation that will be nominally 6-connected
% which hopefully facilitates any subsequent post-processing that seeks to
% merge superpixels.
% - Initial cluster positions are not shifted to point of lowest gradient
% within a 3x3 neighbourhood because this will be rendered irrelevant the
% first time cluster centres are updated.
%
% Reference: R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua and
% S. Susstrunk. "SLIC Superpixels Compared to State-of-the-Art Superpixel
% Methods" PAMI. Vol 34 No 11. November 2012. pp 2274-2281.
%
% See also: SPDBSCAN, MCLEANUPREGIONS, REGIONADJACENCY, DRAWREGIONBOUNDARIES, RGB2LAB
% Copyright (c) 2013 Peter Kovesi
% www.peterkovesi.com
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% Feb 2013
% July 2013 Super pixel attributes returned as a structure array
% Note that most of the computation time is not in the clustering, but rather
% in the region cleanup process.
function [l, Am, Sp, d] = slic(im, k, m, seRadius, colopt, mw, nItr, eim, We)
if ~exist('colopt','var') || isempty(colopt), colopt = 'mean'; end
if ~exist('mw','var') || isempty(mw), mw = 0; end
if ~exist('nItr','var') || isempty(nItr), nItr = 10; end
if exist('eim', 'var'), USEDIST = 0; else, USEDIST = 1; end
MEANCENTRE = 1;
MEDIANCENTRE = 2;
if strcmp(colopt, 'mean')
centre = MEANCENTRE;
elseif strcmp(colopt, 'median')
centre = MEDIANCENTRE;
else
error('Invalid colour centre computation option');
end
[rows, cols, chan] = size(im);
if chan ~= 3
error('Image must be colour');
end
% Convert image to L*a*b* colourspace. This gives us a colourspace that is
% nominally perceptually uniform. This allows us to use the euclidean
% distance between colour coordinates to measure differences between
% colours. Note the image becomes double after conversion. We may want to
% go to signed shorts to save memory.
im = rgb2lab(im);
% Apply median filtering to colour components if mw has been supplied
% and/or non-zero
if mw
if length(mw) == 1
mw(2) = mw(1); % Use same filtering for L and chrominance
end
for n = 1:3
im(:,:,n) = medfilt2(im(:,:,n), [mw(1) mw(1)]);
end
end
% Nominal spacing between grid elements assuming hexagonal grid
S = sqrt(rows*cols / (k * sqrt(3)/2));
% Get nodes per row allowing a half column margin at one end that alternates
% from row to row
nodeCols = round(cols/S - 0.5);
% Given an integer number of nodes per row recompute S
S = cols/(nodeCols + 0.5);
% Get number of rows of nodes allowing 0.5 row margin top and bottom
nodeRows = round(rows/(sqrt(3)/2*S));
vSpacing = rows/nodeRows;
% Recompute k
k = nodeRows * nodeCols;
% Allocate memory and initialise clusters, labels and distances.
C = zeros(6,k); % Cluster centre data 1:3 is mean Lab value,
% 4:5 is row, col of centre, 6 is No of pixels
l = -ones(rows, cols); % Pixel labels.
d = inf(rows, cols); % Pixel distances from cluster centres.
% Initialise clusters on a hexagonal grid
kk = 1;
r = vSpacing/2;
for ri = 1:nodeRows
% Following code alternates the starting column for each row of grid
% points to obtain a hexagonal pattern. Note S and vSpacing are kept
% as doubles to prevent errors accumulating across the grid.
if mod(ri,2), c = S/2; else, c = S; end
for ci = 1:nodeCols
cc = round(c); rr = round(r);
C(1:5, kk) = [squeeze(im(rr,cc,:)); cc; rr];
c = c+S;
kk = kk+1;
end
r = r+vSpacing;
end
% Now perform the clustering. 10 iterations is suggested but I suspect n
% could be as small as 2 or even 1
S = round(S); % We need S to be an integer from now on
for n = 1:nItr
for kk = 1:k % for each cluster
% Get subimage around cluster
rmin = max(C(5,kk)-S, 1); rmax = min(C(5,kk)+S, rows);
cmin = max(C(4,kk)-S, 1); cmax = min(C(4,kk)+S, cols);
subim = im(rmin:rmax, cmin:cmax, :);
assert(numel(subim) > 0)
% Compute distances D between C(:,kk) and subimage
if USEDIST
D = dist(C(:, kk), subim, rmin, cmin, S, m);
else
D = dist2(C(:, kk), subim, rmin, cmin, S, m, eim, We);
end
% If any pixel distance from the cluster centre is less than its
% previous value update its distance and label
subd = d(rmin:rmax, cmin:cmax);
subl = l(rmin:rmax, cmin:cmax);
updateMask = D < subd;
subd(updateMask) = D(updateMask);
subl(updateMask) = kk;
d(rmin:rmax, cmin:cmax) = subd;
l(rmin:rmax, cmin:cmax) = subl;
end
% Update cluster centres with mean values
C(:) = 0;
for r = 1:rows
for c = 1:cols
tmp = [im(r,c,1); im(r,c,2); im(r,c,3); c; r; 1];
C(:, l(r,c)) = C(:, l(r,c)) + tmp;
end
end
% Divide by number of pixels in each superpixel to get mean values
for kk = 1:k
C(1:5,kk) = round(C(1:5,kk)/C(6,kk));
end
% Note the residual error, E, is not calculated because we are using a
% fixed number of iterations
end
% Cleanup small orphaned regions and 'spurs' on each region using
% morphological opening on each labeled region. The cleaned up regions are
% assigned to the nearest cluster. The regions are renumbered and the
% adjacency matrix regenerated. This is needed because the cleanup is
% likely to change the number of labeled regions.
if seRadius
[l, Am] = mcleanupregions(l, seRadius);
else
l = makeregionsdistinct(l);
[l, minLabel, maxLabel] = renumberregions(l);
Am = regionadjacency(l);
end
% Recompute the final superpixel attributes and write information into
% the Sp struct array.
N = length(Am);
Sp = struct('L', cell(1,N), 'a', cell(1,N), 'b', cell(1,N), ...
'stdL', cell(1,N), 'stda', cell(1,N), 'stdb', cell(1,N), ...
'r', cell(1,N), 'c', cell(1,N), 'N', cell(1,N));
[X,Y] = meshgrid(1:cols, 1:rows);
L = im(:,:,1);
A = im(:,:,2);
B = im(:,:,3);
for n = 1:N
mask = l==n;
nm = sum(mask(:));
if centre == MEANCENTRE
Sp(n).L = sum(L(mask))/nm;
Sp(n).a = sum(A(mask))/nm;
Sp(n).b = sum(B(mask))/nm;
elseif centre == MEDIANCENTRE
Sp(n).L = median(L(mask));
Sp(n).a = median(A(mask));
Sp(n).b = median(B(mask));
end
Sp(n).r = sum(Y(mask))/nm;
Sp(n).c = sum(X(mask))/nm;
% Compute standard deviations of the colour components of each super
% pixel. This can be used by code seeking to merge superpixels into
% image segments. Note these are calculated relative to the mean colour
% component irrespective of the centre being calculated from the mean or
% median colour component values.
Sp(n).stdL = std(L(mask));
Sp(n).stda = std(A(mask));
Sp(n).stdb = std(B(mask));
Sp(n).N = nm; % Record number of pixels in superpixel too.
end
%-- dist -------------------------------------------
%
% Usage: D = dist(C, im, r1, c1, S, m)
%
% Arguments: C - Cluster being considered
% im - sub-image surrounding cluster centre
% r1, c1 - row and column of top left corner of sub image within the
% overall image.
% S - grid spacing
% m - weighting factor between colour and spatial differences.
%
% Returns: D - Distance image giving distance of every pixel in the
% subimage from the cluster centre
%
% Distance = sqrt( dc^2 + (ds/S)^2*m^2 )
% where:
% dc = sqrt(dl^2 + da^2 + db^2) % Colour distance
% ds = sqrt(dx^2 + dy^2) % Spatial distance
%
% m is a weighting factor representing the nominal maximum colour distance
% expected so that one can rank colour similarity relative to distance
% similarity. try m in the range [1-40] for L*a*b* space
%
% ?? Might be worth trying the Geometric Mean instead ??
% Distance = sqrt(dc * ds)
% but having a factor 'm' to play with is probably handy
% This code could be more efficient
function D = dist(C, im, r1, c1, S, m)
% Squared spatial distance
% ds is a fixed 'image' we should be able to exploit this
% and use a fixed meshgrid for much of the time somehow...
[rows, cols, chan] = size(im);
[x,y] = meshgrid(c1:(c1+cols-1), r1:(r1+rows-1));
x = x-C(4); % x and y dist from cluster centre
y = y-C(5);
ds2 = x.^2 + y.^2;
% Squared colour difference
for n = 1:3
im(:,:,n) = (im(:,:,n)-C(n)).^2;
end
dc2 = sum(im,3);
D = sqrt(dc2 + ds2/S^2*m^2);
%--- dist2 ------------------------------------------
%
% Usage: D = dist2(C, im, r1, c1, S, m, eim)
%
% Arguments: C - Cluster being considered
% im - sub-image surrounding cluster centre
% r1, c1 - row and column of top left corner of sub image within the
% overall image.
% S - grid spacing
% m - weighting factor between colour and spatial differences.
% eim - Edge strength sub-image corresponding to im
%
% Returns: D - Distance image giving distance of every pixel in the
% subimage from the cluster centre
%
% Distance = sqrt( dc^2 + (ds/S)^2*m^2 )
% where:
% dc = sqrt(dl^2 + da^2 + db^2) % Colour distance
% ds = sqrt(dx^2 + dy^2) % Spatial distance
%
% m is a weighting factor representing the nominal maximum colour distance
% expected so that one can rank colour similarity relative to distance
% similarity. try m in the range [1-40] for L*a*b* space
%
function D = dist2(C, im, r1, c1, S, m, eim, We)
% Squared spatial distance
% ds is a fixed 'image' we should be able to exploit this
% and use a fixed meshgrid for much of the time somehow...
[rows, cols, chan] = size(im);
[x,y] = meshgrid(c1:(c1+cols-1), r1:(r1+rows-1));
x = x-C(4);
y = y-C(5);
ds2 = x.^2 + y.^2;
% Squared colour difference
for n = 1:3
im(:,:,n) = (im(:,:,n)-C(n)).^2;
end
dc2 = sum(im,3);
% Combine colour and spatial distance measure
D = sqrt(dc2 + ds2/S^2*m^2);
% for every pixel in the subimage call improfile to the cluster centre
% and use the largest value as the 'edge distance'
rCentre = C(5)-r1; % Cluster centre coords relative to this sub-image
cCentre = C(4)-c1;
de = zeros(rows,cols);
for r = 1:rows
for c = 1:cols
v = improfile(eim,[c cCentre], [r rCentre]);
de(r,c) = max(v);
end
end
% Combine edge distance with weight, We with total Distance.
D = D + We * de;