DIP数字图像处理（二）超像素 SLIC简单线性迭代聚类 matlab实现（不调库）

SLIC观测常用五维向量$$\text{z}=(r,g,b,x,y)$$
令$n_{sp}$是所需要的超像素数，$n_{tp}$是图像中的像素总数。初始的超像素中心${\text{m}}_i=[r_i,g_i,b_i,x_i,y_i{]}^T$，$i=1,2,\cdots n_{sp}$是对图像以间距为$s$个单位的规则网格取样得到的。选择$s=[n_{tp}/n_{sp}{]}^{1/2}$，初始聚类中心设为每个中心周围$3\times 3$邻域的最小梯度位置。

算法步骤

1. 初始化算法：以$s$对图像取样，计算初始超像素聚类中心$${\text{m}}_i=[r_i,g_i,b_i,x_i,y_i{]}^T$$将聚类中心移到$3\times 3$邻域的最小梯度位置。对于图像中每个像素位置$p$，设置$L(p)=-1$和距离$d(p)=inf$.
   （未完待续）

代码实现（不调库）

clear;clc;close all;
img = imread("./images/sea_house.jpg");
imshow(img,[]);
[m,n,~]=size(img);

n_tp = (m-2)*(n-2);
n_sp = 1000;
s = round((n_tp/n_sp)^0.5);
M = zeros(5,10);

% 初始聚类中心
index = 1;
for i = 2:s:m-1
    for j = 2:s:n-1
        M(:,index)=[img(i,j,1),img(i,j,2),img(i,j,3),i,j]';
        M(4,index)=i;
        M(5,index)=j;
        index=index+1;
    end
end
label=zeros(m,n);
label=label-1;
d=zeros(m,n);
d=d+999999;

[~, Mcol] = size(M);
temp=zeros(m,n);
for i=1:Mcol
    temp(M(4,i),M(5,i))=1;
end
imshow(temp,[]);

% 计算梯度
img_d=double(img)/255;
r=img_d(:,:,1);
g=img_d(:,:,2);
b=img_d(:,:,3);
Y=0.299*r+0.587*g+0.114*b;

f1=fspecial('sobel');
f2=f1';
gx=imfilter(Y,f1);
gy=imfilter(Y,f2);
G=sqrt(gx.^2+gy.^2);

% 分配聚类中心
[~, Mcol] = size(M);
for index = 1:Mcol
    x = M(4,index);
    y = M(5,index);
    grad_current = img_d(x,y);
    x_current = x;
    y_current = y;
    for x_new = x-1:x+1
        for y_new = y-1:y+1
            if img_d(x_new,y_new)<grad_current
                grad_current = img_d(x_new,y_new);
                x_current = x_new;
                y_current = y_new;
            end
        end
    end
    M(:,index) = [img(x_current,y_current,1),img(x_current,y_current,2),img(x_current,y_current,3),x_current,y_current]';
    M(4,index) = x_current;
    M(5,index) = y_current;
end

% 将样本分配给聚类中心
%flag=0;
for iter = 1:100
vec = zeros(5,1);
for index = 1:Mcol
    x = M(4,index);
    y = M(5,index);
    for i = x-s:x+s-1
        for j = y-s:y+s-1
            if i<1||j<1||i>m||j>n
                continue;
            end
            vec(1,1) = img(i,j,1);
            vec(2,1) = img(i,j,2);
            vec(3,1) = img(i,j,3);
            vec(4,1) = i;
            vec(5,1) = j;
%             if vec(4,1)>255
%                 flag=1;
%             end
            distance = norm(M(:,index)-double(vec));
            if distance < d(i,j)
                d(i,j)=distance;
                label(i,j)=index;
            end
        end
    end
end

% update the center
M_new = zeros(5,Mcol);
counter = zeros(1,Mcol);
for i = 1:m
    for j = 1:n
        l = label(i,j);
        vec_n = zeros(5,1);
        vec_n(1,1) = img(i,j,1);
        vec_n(2,1) = img(i,j,2);
        vec_n(3,1) = img(i,j,3);
        vec_n(4,1) = i;
        vec_n(5,1) = j;
        M_new(:,l) = M_new(:,l) + vec_n;
        counter(1,l) = counter(1,l)+1;
    end
end

% convergence check
for i = 1:Mcol
    M_new(:,i) = M_new(:,i)./counter(1,i);
end
diff = norm(M_new-M);
if diff<10
    break
end

M = round(M_new);
end

% average
img_final = zeros(m,n,3);
for i = 1:m
    for j = 1:n
        l = label(i,j);
        img_final(i,j,1) = M(1,l);
        img_final(i,j,2) = M(2,l);
        img_final(i,j,3) = M(3,l);
    end
end
figure
imshow(uint8(img_final),[]);
title("center = 100");