Spectral Matting的python实现

Spectral Matting是一个很好的矩阵分析、图像处理例子，是一个很好的学习例子。
在Github上，有Spectral Matting的MATLAB的实现，却没找到Python的实现，于是我参考了以下两个代码源：
[1]https://github.com/yaksoy/SemanticSoftSegmentation
[2]https://github.com/MarcoForte/closed-form-matting
用python实现了一次，这次实践过程收益良多。

[1]是语义软分割（Yagiz Aksoy, Tae-Hyun Oh, Sylvain Paris, Marc Pollefeys and Wojciech Matusik, “Semantic Soft Segmentation”, ACM Transactions on Graphics (Proc. SIGGRAPH), 2018）的MATLAB实现，它的实现参考了Spectral Matting 原MATLAB的实现方式(http://www.vision.huji.ac.il/SpectralMatting/)，并在其上改动了一些部分，让代码更为清晰，为此我选了这个版本作为python翻译的参考。
[2]是A closed-form solution to natural image matting的python实现，Spectral Matting其实是closed-form matting的后续，它们在Laplacian矩阵生成的方式上，是完成相同的。
接下来，我将分以下几部分来实现Spectral Matting：

• Laplician矩阵的构建
• 特征值与特征向量选择
• 优化——寻找合适的alpha
• alpha的初始化
• 优化迭代的实现
• Components的组合

---

1、Laplician矩阵的构建

根据[1]，图像的Laplicain矩阵中元素为：
$$A_q(i,j)=\{\begin{array}{cc}{\delta }_{ij}-\frac{1}{|w_q|}\left(1+(I_i-u_q{)}^T({\Sigma }_q+\frac{ϵ}{|w_q|}{I}_{3\times 3})(I_j-u_q)\right)& i,j\in w_q\\ 0& \text{otherwise}\end{array}(1)$$
其中，$w_q$是$3\times 3$窗体，$|w_q|$是窗体的pixels的数量，${\Sigma }_q$是该窗体的$3\times 3$协方差，$I_i,I_j$是窗体中像素i和像素j的颜色值。
这部分的实现代码如下：

def _rolling_block(A, block=(3, 3)):
    """Applies sliding window to given matrix."""
    shape = (A.shape[0] - block[0] + 1, A.shape[1] - block[1] + 1) + block
    strides = (A.strides[0], A.strides[1]) + A.strides
    return as_strided(A, shape=shape, strides=strides)
    
def mattingAffinity(img, eps=10**(-7), win_rad=1):
    """Computes Matting Laplacian for a given image.
    Args:
        img: 3-dim numpy matrix with input image
        eps: regularization parameter controlling alpha smoothness
            from Eq. 12 of the original paper. Defaults to 1e-7.
        win_rad: radius of window used to build Matting Laplacian (i.e.
            radius of omega_k in Eq. 12).
    Returns: sparse matrix holding Matting Laplacian.
    """
    win_size = (win_rad * 2 + 1) ** 2
    h, w, d = img.shape
    # Number of window centre indices in h, w axes
    c_h, c_w = h - 2 * win_rad, w - 2 * win_rad
    win_diam = win_rad * 2 + 1

    indsM = np.arange(h * w).reshape((h, w))
    ravelImg = img.reshape(h * w, d)
    win_inds = _rolling_block(indsM, block=(win_diam, win_diam))

    #win_inds = win_inds.reshape(c_h, c_w, win_size)
    
    win_inds = win_inds.reshape(-1, win_size)
   
    winI = ravelImg[win_inds]

    win_mu = np.mean(winI, axis=1, keepdims=True)
    win_var = np.einsum('...ji,...jk ->...ik', winI, winI) / win_size - np.einsum('...ji,...jk ->...ik', win_mu, win_mu)

    inv = np.linalg.inv(win_var + (eps/win_size)*np.eye(3))

    X = np.einsum('...ij,...jk->...ik', winI - win_mu, inv)
    X = np.einsum('...ij,...kj->...ik', X, winI - win_mu)
    vals = (1.0/win_size)*(1 + X)

    nz_indsCol = np.tile(win_inds, win_size).ravel()
    nz_indsRow = np.repeat(win_inds, win_size).ravel()
    nz_indsVal = vals.ravel()
    A = scipy.sparse.coo_matrix((nz_indsVal, (nz_indsRow, nz_indsCol)), shape=(h*w, h*w))
    return A

def affinityMatrixToLaplacian(aff):
    degree = aff.sum(axis=1)
    degree = np.asarray(degree)
    degree = degree.ravel()
    Lap = scipy.sparse.diags(degree)-aff
    return Lap

其中_rolling_block可将图像按窗口要求（3*3）切成小块，并返回窗体像素所在矩阵的index，用于窗口的像素提取。
公式(1)的具体实现时以下代码片段：

win_mu = np.mean(winI, axis=1, keepdims=True)
    win_var = np.einsum('...ji,...jk ->...ik', winI, winI) / win_size - np.einsum('...ji,...jk ->...ik', win_mu, win_mu)

    inv = np.linalg.inv(win_var + (eps/win_size)*np.eye(3))

    X = np.einsum('...ij,...jk->...ik', winI - win_mu, inv)
    X = np.einsum('...ij,...kj->...ik', X, winI - win_mu)
    vals = (1.0/win_size)*(1 + X)

此处用到np.einsum()函数，它可以完成复杂的张量计算。在得到各窗体计算结果后，由scipy.sparse.coo_matrix()构造稀疏矩阵L。
将前面几个函数串起来，求图像的laplacian矩阵由以下代码片段实现，其中图像来自：
https://github.com/yaksoy/SemanticSoftSegmentation/blob/master/docia.png

image = Image.open('docia.png')
img_s = image.resize((160,80))
im_np = np.asarray(img_s)
w = im_np.shape[0]
img = im_np[:,:w,:]/255.0

A = mattingAffinity(img)
L = affinityMatrixToLaplacian(A)

图1 原图
图像读出来后转换为numpy的array类型，将它作为mattingAffinity()的输入，计算亲和矩阵A，和Laplician矩阵L。

2、特征值与特征向量选择

[1]的一个关键思路是：在图像分割后，能得到若干个components（比如是10个），它们可由最小特征矢量（比如是50个）的线性组合得到，即：
$${\alpha }^k=\mathbf{E}{\mathbf{y}}^k(2)$$
其中，$\mathbf{E}$ 是特征矢量（比如是，最小的50个特征矢量）构成的矩阵，为对应我们的上述代码，其形状为$6400\times 50$，6400是图像像素的总数（$80\times 80$），50是最小特征矢量的数量。${\mathbf{y}}^k$ 是第k个component对应的线性组合，其形状为$50\times 1$，而 α k \alpha^k