SuperPixel

Gonzalez R. C. and Woods R. E. Digital Image Processing (Forth Edition).

单个像素的意义其实很小, 于是有了superpixel的概念, 即一簇pixels的集合(且这堆pixels共用一个值), 这会导致图片有非常有趣的艺术风格(下图便是取不同的superpixel大小形成的效果, 有种抽象画的感觉?):

经过superpixel的预处理后, 图片可以变得更加容易提取edge, region, 毕竟superpixel已经率先提取过一次了.

SLIC Superpixel algorithm

SLIC (simple linear iterative clustering) 算法是基于k-means的一种聚类算法.

Given: 需要superpixels的个数$n_{sp}$; 图片$f(x,y)=(r,g,b),x=1,2,\cdots M,y=1,2,\cdots ,N$;

1. 根据图片以及其位置信息生成数据:
   $$\mathbf{z}=[r,g,b,x,y{]}^T,$$
   其中$r,g,b$是颜色编码, $x,y$是位置信息.

2. 令$n_{tp}=MN$表示pixels的个数, 并计算网格大小:
   $$s=[n_{tp}/n_{sp}{]}^{1/2}.$$

3. 将图片均匀分割为大小$s$的网格, 初始化superpixels的中心:
   $${\mathbf{m}}_i=[r_i,g_i,b_i,x_i,y_i{]}^T,i=1,2,\cdots ,n_{sp},$$
   为网格的中心. 或者, 为了防止噪声的影响, 选择中心$3\times 3$领域内梯度最小的点.

4. 将图片的每个pixel的类别标记为$L(p)=-1$, 距离$d(p)=\infty$;

5. 重复下列步骤直到收敛:

   1. 对于每个像素点$p$, 计算其与$2s\times 2s$邻域内的中心点${\mathbf{m}}_i$之间的距离$D_i(p)$, 倘若$D_i(p)<d(p)$:
      $$d(p)=D_i,L(p)=i.$$

   2. 令$C_i$表示$L(p)=i$的像素点的集合, 更新superpixels的中心:
      $${\mathbf{m}}_i=\frac{1}{|C_i|}\sum _{\mathbf{z}\in C_i}\mathbf{z},i=1,2,\cdots ,n_{sp}.$$

6. 将以${\mathbf{m}}_i$为中心的区域中的点的(r, g, b)设定为与${\mathbf{m}}_i$一致.

距离函数的选择

倘若$D$采用的是和普通K-means一样的$\Vert \cdot {\Vert }_2$显然是不合适的, 因为$(r,g,b)$和$(x,y)$显然不是一个尺度的. 故采用如下的距离函数:
$$\begin{gathered} D=[(\frac{d_c}{d_{cm}}{)}^2+(\frac{d_s}{d_{sm}}{)}^2{]}^{1/2}, \\ {d}_c=[(r_j-r_i{)}^2+(g_j-g_i{)}^2+(b_j-b_i{)}^2{]}^{1/2}, \\ {d}_s=[(x_j-x_i{)}^2+(y_j-y_i{)}^2{]}^{1/2}, \end{gathered}$$
其中$d_{cm},d_{sm}$分别是$d_c,d_s$可能取到的最大值, 相当于标准化了.

代码

skimage.segmentation.slic

import numpy as np


def _generate_data(img):
    img = img.astype(np.float64)
    if len(img.shape) == 2:
        img = img[..., None]
    M, N = img.shape[0], img.shape[1]
    loc = np.stack(np.meshgrid(range(M), range(N), indexing='ij'), axis=-1)
    classes = -np.ones((M, N))
    distances = np.ones((M, N)) * np.float('inf')
    data = np.concatenate((img, loc), axis=-1)
    return data, classes, distances

def _generate_means(data, size: int):
    M, N = data.shape[0], data.shape[1]
    x_splits = np.arange(0, M + size, size)
    y_splits = np.arange(0, N + size, size)
    means = []
    for i in range(len(x_splits) - 1):
        for j in range(len(y_splits) - 1):
            r1, r2 = x_splits[i:i+2]
            c1, c2 = y_splits[j:j+2]
            region = data[r1:r2, c1:c2]
            means.append(region.mean(axis=(0, 1)))
    return np.array(means)


def _unit_step(data, means, classes, distances, size, dis_fn):
    M, N = data.shape[0], data.shape[1]
    size = 2 * size
    for i, m in enumerate(means):
        # ..., x, y
        x, y = np.round(m[-2:])
        x, y = int(x), int(y)
        xl, xr = max(0, x - size), min(x + size, M)
        yb, yt = max(0, y - size), min(y + size, N)
        p = data[xl:xr, yb:yt]
        _dis = dis_fn(p, m)
        indices = _dis < distances[xl:xr, yb:yt]
        distances[xl:xr, yb:yt][indices] = _dis[indices]
        classes[xl:xr, yb:yt][indices] = i

    # update
    for i in range(len(means)):
        x_indices, y_indices = np.where(classes == i)
        if len(x_indices) == 0:
            continue
        means[i] = data[x_indices, y_indices].mean(axis=0)

def slic(img, size, max_iters=10, compactness=10):
    data, classes, distances = _generate_data(img)
    means = _generate_means(data, size)
    dsm = size
    dcm = (img.max(axis=(0, 1)) - img.min(axis=(0, 1))) * compactness
    dsc = np.concatenate((dcm, [dsm] * 2))
    def dis_func(p, m):
        _dis = ((p - m) / dsc) ** 2
        return _dis.sum(axis=-1)
    for _ in range(max_iters):
        _unit_step(data, means, classes, distances, size, dis_func)
    new_img = np.zeros_like(img, dtype=np.float)
    for i, m in enumerate(means):
        x_indices, y_indices = np.where(classes == i)
        if len(x_indices) == 0:
            continue
        new_img[x_indices, y_indices] = m[:-2]
    return new_img.astype(img.dtype)
    

from skimage import io, segmentation, filters
from freeplot.base import FreePlot

    
img = io.imread(r"Lenna.png")    

ours = slic(img, size=50, compactness=0.5)

def mask2img(mask, img):
    new_img = img.astype(np.float)
    masks = np.unique(mask)
    for m in masks:
        x, y = np.where(mask == m)
        mcolor = new_img[x, y].mean(axis=0)
        new_img[x, y] = mcolor
    return new_img.astype(img.dtype)


mask = segmentation.slic(img)
yours = mask2img(mask, img)

fp = FreePlot((1, 3), (10.3, 5), titles=('Lenna', 'ours', 'skimage.segmentation.slic'))
fp.imageplot(img, index=(0, 0))
fp.imageplot(ours, index=(0, 1))
fp.imageplot(yours, index=(0, 2))
fp.set_title()
fp.show()

skimage上实现的代码还有强制连通性, 我想这个是为什么它看起来这么流畅的原因. Compactness 越大, 聚类越倾向于空间信息, 所以越容易出现块状结构.