一、簡介

1 概念
超畫素由一系列位置相鄰且顏色、亮度、紋理等特徵相似的畫素點組成的小區域。這些小區域大多保留了進一步進行影象分割的有效資訊，且一般不會破壞影象中物體的邊界資訊。
超畫素是吧一幅畫素級（pixel-level）的圖，劃分成區域級（district-level）的圖，是對基本資訊元素進行的抽象。

（a）是原始影象，（b）是基於人類視角的分割圖（groundtruth），（c）是超畫素分割的影象，（d）是基於（c）進行分割的影象。
超畫素最大的功能之一是作為影象處理其他演算法的預處理，在不犧牲太大精確度的情況下降維。
超畫素最直觀的解釋是把一些具有相似特性的畫素“聚合”起來，形成一個更具有代表性的大“元素”。而這個新元素，將作為其他影象處理演算法的基本單位。這樣可以降低維度，剔除一些異常畫素點。
理論上，任何影象分割演算法的過度分割（over-segmentation）即可生成超畫素。
影象分割中的超畫素是指具有相似紋理、顏色、亮度等特徵的相鄰相似構成的具有一定意義的不規則的畫素塊。它利用畫素之間特徵的相似性將畫素分組，用少量的超畫素代替大量的畫素來表達影象特徵，很大程度上降低了影象處理的複雜度，所以通常作為分割演算法的預處理步驟。

2 超畫素判別條件
（1）Undersegmentation Error

（2）如上圖所示，白色是影象中的一個物體，紅線是一個個超畫素的輪廓，而粉紅色的區域就是undersegmentation區域，這部分割槽域越大越不好。
Boundary Recall

如上圖所示，黑色虛線及實線是影象中物體的輪廓，紅線是超畫素的邊界。一個好的超畫素演算法，應該覆蓋影象中物體的輪廓。在給予一定緩衝（粉紅色區域）的情況下，超畫素邊緣可以覆蓋影象物體邊緣（黑色實線）越多越好。
（3）Compactness score
衡量超畫素是否“緊實”。

3 超畫素初始化的方法
種子畫素初始化
SLIC利用了簡單的聚類（貪婪）演算法，初始時，每一個聚類的中心被平均的分佈在影象中，而超畫素的個數，可以基本由這些中心點來決定。每一步迭代，種子畫素合併周圍的畫素，形成超畫素。

矩形區域初始化
SEEDS的初始化是把影象平均分割成很多矩形，初始超畫素即為這些矩形。每一步迭代，超畫素的邊緣不斷變化，直到匯合。

4 超畫素演算法

5 SLIC演算法
SLIC（simple linear iterative clustering），即簡單的線性迭代聚類。它是2010年提出的一種思想簡單、實現方便的演算法，將彩色影象轉換為CIELAB顏色空間和XY座標下的5維特徵向量，然後對5維特徵向量構造距離度量標準，對影象畫素進行區域性聚類的過程。SLIC演算法能生成緊湊近似均勻的超畫素，在運算速度，物體輪廓保持、超畫素形狀方面具有較高的綜合評價，比較符合人們期望的分割效果。

5.1 SLIC優點：

生成的超畫素如同細胞一般緊湊整齊，鄰域特徵比較容易表達。這樣基於畫素方法可以比較容易的改造為基於超畫素的方法。
不僅可以分割彩色影象，也可以相容分割灰度圖。
需要設定的引數非常少，預設情況下只需要設定一個預分割的超畫素的數量。
相比其他的超畫素的分割方法，SLIC在執行速度、生成超畫素的緊湊度、輪廓保持方面都比較理想。

5.2 演算法步驟：

初始化種子點（聚類中心）：按照設定的超畫素的個數，在影象內均勻的分配種子點。假設影象總共有N個畫素點，預分割為K個相同尺寸的超畫素，那麼每個超畫素的大小為N/K，則相鄰種子點的距離（步長）近似為S=sqrt(N/K)。
在種子點的n*n領域內重新選擇種子點（一般取n=3）：計算該領域內所有畫素點的梯度值，將種子點移到該領域內梯度最小的地方。避免種子點落在梯度較大的輪廓邊界上，以免影響後續聚類效果。
在每個種子點周圍的領域內為每個畫素點分配類別標籤（即屬於哪個聚類中心）：SLIC的搜尋範圍是2S2S，期望的超畫素尺寸為SS，這樣可以加速演算法收斂。
距離度量：包括顏色距離和空間距離。對每個搜尋到的畫素點，分別計算它和該種子點的距離。

二、原始碼

% 能夠確定超畫素位置的核心程式是  EnforceLabelConnectivity
% 11-123頁碼 
clc
clear
close all;
tic
img = imread('bee.jpg');
imshow(img)
title('original')
%設定超畫素個數
K = 500;
%設定超畫素緊湊係數
m_compactness = 100;

%%
img = DeSample(img,2);
img_size = size(img);
%轉換到LAB色彩空間
cform = makecform('srgb2lab');       %rgb空間轉換成lab空間 matlab自帶的用法,Create color transformation structure
img_Lab = applycform(img, cform);    %rgb轉換成lab空間
figure;
imshow(img_Lab)
title('img_lab')

%%
%得到超畫素的LABXY種子點資訊     
img_sz = img_size(1)*img_size(2);
superpixel_sz = img_sz/K;  % 每個超畫素的畫素點數
STEP = uint32(sqrt(superpixel_sz)); % 開方的邊長
xstrips = uint32(img_size(2)/STEP);  % x方向 的超畫素個數
ystrips = uint32(img_size(1)/STEP);  % y方向 的超畫素個數
xstrips_adderr = double(img_size(2))/double(xstrips);  
ystrips_adderr = double(img_size(1))/double(ystrips);
numseeds = xstrips*ystrips;   % 實際的超畫素個數
%種子點xy資訊初始值為晶格中心亞畫素座標
%種子點Lab顏色資訊為對應點最接近畫素點的顏色通道值
kseedsx = zeros(numseeds, 1);
kseedsy = zeros(numseeds, 1);
kseedsl = zeros(numseeds, 1);
kseedsa = zeros(numseeds, 1);
kseedsb = zeros(numseeds, 1);

n = 1;
for y = 1: ystrips   % 第y個超畫素
    for x = 1: xstrips   % 第 x 個超畫素
        kseedsx(n, 1) = (double(x)-0.5)*xstrips_adderr; % 第x個種子點中心座標，非準確描述
        kseedsy(n, 1) = (double(y)-0.5)*ystrips_adderr; % 第y個種子點中心座標，非準確描述
        % 種子點中心對應LAB圖上位置的 三通道值
        kseedsl(n, 1) = img_Lab(fix(kseedsy(n, 1)), fix(kseedsx(n, 1)), 1);  % fix 417.1296 變417
        kseedsa(n, 1) = img_Lab(fix(kseedsy(n, 1)), fix(kseedsx(n, 1)), 2); 
        kseedsb(n, 1) = img_Lab(fix(kseedsy(n, 1)), fix(kseedsx(n, 1)), 3);
        n = n+1;
    end
end

n = 1;
%根據種子點計算超畫素分割槽
klabels = PerformSuperpixelSLIC(img_Lab, kseedsl, kseedsa, kseedsb, kseedsx, kseedsy, STEP, m_compactness);
%合併小的分割槽
[supmtrx,supmtry,nlabels] = EnforceLabelC(img_Lab, klabels, K);
% 這裡的supmtrx,supmtry的每列分別是對應標籤區域的全部x座標和y座標
function klabels = PerformSuperpixelSLIC(img_Lab, kseedsl, kseedsa, kseedsb, kseedsx, kseedsy, STEP, compactness)

[m_height, m_width, m_channel] = size(img_Lab);
[numseeds xxxxx]= size(kseedsl);
img_Lab = double(img_Lab);
%畫素標籤格式為(x, y) (行, 列)
klabels = zeros(m_height, m_width);
%聚類尺寸
clustersize = zeros(numseeds,1);
inv = zeros(numseeds,1);
sigmal = zeros(numseeds,1);
sigmaa = zeros(numseeds,1);
sigmab = zeros(numseeds,1);
sigmax = zeros(numseeds,1);
sigmay = zeros(numseeds,1);
invwt = 1/( (double(STEP)/double(compactness)) *(double(STEP)/double(compactness)) );
%invwt = double(compactness)/double(STEP);
distvec = 100000*ones(m_height, m_width);
numk = numseeds;
for itr = 1: 10   %迭代次數
    sigmal = zeros(numseeds, 1);
    sigmaa = zeros(numseeds, 1);
    sigmab = zeros(numseeds, 1);
    sigmax = zeros(numseeds, 1);
    sigmay = zeros(numseeds, 1);
    clustersize = zeros(numseeds, 1);
    inv = zeros(numseeds, 1);
    distvec = double(100000*ones(m_height, m_width));
    %根據當前種子點資訊計算每一個畫素的歸屬
    for n = 1: numk
        y1 = max(1, kseedsy(n, 1)-STEP);
        y2 = min(m_height, kseedsy(n, 1)+STEP);
        x1 = max(1, kseedsx(n, 1)-STEP);
        x2 = min(m_width, kseedsx(n, 1)+STEP);
        %按畫素計算距離
        for y = y1: y2
            for x = x1: x2
                %dist_lab = abs(img_Lab(y, x, 1)-kseedsl(n))+abs(img_Lab(y, x, 2)-kseedsa(n))+abs(img_Lab(y, x, 3)-kseedsb(n));
                % lab圖 點到種子點 定義距離差，判斷相似度
                dist_lab = (img_Lab(y, x, 1)-kseedsl(n, 1))^2+(img_Lab(y, x, 2)-kseedsa(n, 1))^2+(img_Lab(y, x, 3)-kseedsb(n, 1))^2;
                % 改成平方啊 ！！！  @@@  ！！！  @@@   ！！！
                dist_xy = (double(y)-kseedsy(n, 1))*(double(y)-kseedsy(n, 1)) + (double(x)-kseedsx(n, 1))*(double(x)-kseedsx(n, 1));
                %dist_xy = abs(y-kseedsy(n)) + abs(x-kseedsx(n));
                
                %距離 = lab色彩空間距離 + 空間距離權重×空間距離
                dist = dist_lab + dist_xy*invwt;
                %在周圍最多四個種子點中找到最相似的 標記後存入klabels
                %m = (y-1)*m_width+x;
                if (dist<distvec(y, x))
                    distvec(y, x) = dist;  % 不斷變小
                    klabels(y, x) = n;  % n是標籤，也就是值畫素屬於哪個種子點
                end
            end
        end
    end
    %完成一遍分類後，重新計算種子點位置 使其向梯度最小地方移動
    ind = 1;
    for r = 1: m_height
        for c = 1: m_width
            sigmal(klabels(r, c),1) = sigmal(klabels(r, c),1)+img_Lab(r, c, 1);  % 畫素塊內所有的通道值相加
            sigmaa(klabels(r, c),1) = sigmaa(klabels(r, c),1)+img_Lab(r, c, 2);
            sigmab(klabels(r, c),1) = sigmab(klabels(r, c),1)+img_Lab(r, c, 3);
            sigmax(klabels(r, c),1) = sigmax(klabels(r, c),1)+c;    % 畫素塊內所有的橫座標相加
            sigmay(klabels(r, c),1) = sigmay(klabels(r, c),1)+r;     % 畫素塊內所有的縱座標相加
            clustersize(klabels(r, c),1) = clustersize(klabels(r, c),1)+1;   % 畫素塊內所有個數相加
        end
    end
    for m = 1: numseeds  % 第m個種子點
        if (clustersize(m, 1)<=0)
            clustersize(m, 1) = 1;
        end
        inv(m, 1) = 1/clustersize(m, 1);
    end
    function [supmtrx,supmtry,nlabels] = EnforceLabelC(img_Lab, labels, K)

dx = [-1, 0, 1, 0];  %四鄰域
dy = [0, -1, 0, 1];
[m_height, m_width, m_channel] = size(img_Lab);
[M, N] = size(labels);
numlabel = max(max(labels));
SUPSZ = (m_height*m_width)/K;   %標準區域面積
nlabels = (-1)*ones(M, N);

label = 1;
adjlabel = 1;
xvec = zeros(m_height*m_width, 1);
yvec = zeros(m_height*m_width, 1);
supmtrx = zeros(2*floor(SUPSZ), numlabel);
supmtry = zeros(2*floor(SUPSZ), numlabel);
m = 1;
n = 1;

for j = 1: m_height
    for k = 1: m_width
        %逐點尋找未標記的區域 小於0 才執行
        if (0>nlabels(m, n))
            %從第一個未標記的(m,n)起，確定一個新區域，用label標記該區域的起點，用蝶形前進
            nlabels(m, n) = label;
            %開始一個新的分割 記錄起點座標
            xvec(1, 1) = k;
            yvec(1, 1) = j;
            supmtrx(1, label) = k;
            supmtry(1, label) = j;
            %如果起點與某個已知區域相連 用adjlabel記錄該區域編號 如果當前區域過小則與相鄰區域合併
            for i = 1: 4
                x = xvec(1, 1)+dx(1, i);
                y = yvec(1, 1)+dy(1, i);
                if (x>0 && x<=m_width && y>0 && y<=m_height)
                    if (nlabels(y, x)>0)
                        adjlabel = nlabels(y, x);  % 一般是左臨或上鄰的標籤
                    end
                end
            end

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423