% 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, espiecially 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
% Centre for Exploration Targeting
% School of Earth and Environment
% The University of Western Australia
% peter.kovesi at uwa edu au
%
% 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);
