匹配追踪算法进行图像重建

匹配追踪的过程已经在匹配追踪算法（MP）简介中进行了简单介绍，下面是使用Python进行图像重建的实践。

MP算法Python版

MP算法原理：

算法假定输入信号与字典库中的原子在结构上具有一定的相关性，这种相关性通过信号与原子库中原子的内积表示，即内积越大，表示信号与字典库中的这个原子的相关性越大，因此可以使用这个原子来近似表示这个信号。当然这种表示会有误差，将表示误差称为信号残差，用原信号减去这个原子，得到残差，再通过计算相关性的方式从字典库中选出一个原子表示这个残差。迭代进行上述步骤，随着迭代次数的增加，信号残差将越来越小，当满足停止条件时终止迭代，得到一组原子，及残差，将这组原子进行线性组合就能重构输入信号。

MP算法的执行步骤如下：

输入：字典矩阵A，信号向量y，稀疏度k.

输出：x的k稀疏逼近x^.

初始化：生成字典矩阵A（这里使用离散余弦变换基DCT），残差r0=y，索引集Λ0=∅，t=1.

循环执行步骤1-5：

1. 找出残差r和字典矩阵的列Ai积中最大值所对应的值p及脚标λ，即pt=maxi=1,⋯,N|<rt−1,Ai>|.
2. 更新索引集Λt=Λt−1∪{λt}，记录找到的字典矩阵中的重建原子集合At=[At−1,Aλt].
3. 更新稀疏向量x^t=x^t∪{pt}.
4. 更新残差rt=y−Atx^t，t=t+1.
5. 判断是否满足t>k，若满足，则迭代停止；若不满足，则继续执行步骤1.

Python代码实现（针对二维图像）：

import numpy as npdef bmp(mtx, codebook, threshold):    """    :param mtx: 原始图像（mxn）    :param codebook: 字典（mxk）    :param threshold: 非零元素个数的最大值    :return: 稀疏编码系数   3 """    n = mtx.shape[1] if len(mtx.shape) > 1 else 1  # 原始图像mtx中向量的个数    k = codebook.shape[1]  # 字典dictionary中向量的个数    result = np.zeros((k, n))  # 系数矩阵result中行数等于dictionary中向量的个数，列数等于mtx中向量的个数    for i in range(n):        indices = []  # 记录选中字典中原子的位置        coefficients = []  # 存储系数向量        residual = mtx[:, i]        for j in range(threshold):            projection = np.dot(codebook.T, residual)            # 获取内积向量中元素绝对值的最大值            max_value = projection.max()            if abs(projection.min()) >= abs(projection.max()):                max_value = projection.min()            pos = np.where(projection == max_value)[0]            indices.append(pos.tolist()[0])  # 只存储在字典中的列（因为计算过程中对codebook进行了转置，所以这里取第一个元素）            coefficients.append(max_value)            residual = mtx[:, i] - np.dot(codebook[:, indices[0: j + 1]], np.array(coefficients))            if (residual ** 2).sum() < 1e-6:                break        for t, s in zip(indices, coefficients):            result[t][i] = s    return result

基于MP的图像重建

对于较大的图像，进行分块处理，使用im2col和col2im函数进行图像的分块和分块后的重建（参考：Python中如何实现im2col和col2im函数）。

这样字典矩阵的行数就仅仅和分块矩阵的大小有关，和原始图像的大小没有关系了。我们可以使用规模较小的字典矩阵表征较大的图像。

Python代码实现：

import numpy as npfrom scipy import fftpackimport mathimport mahotas as mhimport matplotlib.pyplot as pltimport mp.mpalgdef dct2(mtx):    return fftpack.dct(fftpack.dct(mtx.T, norm='ortho').T, norm='ortho')def idct2(mtx):    return fftpack.idct(fftpack.idct(mtx.T, norm='ortho').T, norm='ortho')def dctmtx(n):    basis = np.zeros((n, n))    for i in range(n):        c = math.sqrt(2 / n) if i != 0 else math.sqrt(1 / n)        for j in range(n):            basis[i, j] = c * math.cos((j + 0.5) * math.pi * i / n)    return basisdef im2col(mtx, block_size):    mtx_shape = mtx.shape    sx = mtx_shape[0] - block_size[0] + 1    sy = mtx_shape[1] - block_size[1] + 1    # 如果设A为m×n的，对于[p q]的块划分，最后矩阵的行数为p×q，列数为(m−p+1)×(n−q+1)。    result = np.empty((block_size[0] * block_size[1], sx * sy))    # 沿着行移动，所以先保持列（i）不动，沿着行（j）走    for i in range(sy):        for j in range(sx):            result[:, i * sx + j] = mtx[j:j + block_size[0], i:i + block_size[1]].ravel(order='F')    return resultdef col2im(mtx, image_size, block_size):    p, q = block_size    sx = image_size[0] - p + 1    sy = image_size[1] - q + 1    result = np.zeros(image_size)    weight = np.zeros(image_size)  # weight记录每个单元格的数字重复加了多少遍    col = 0    # 沿着行移动，所以先保持列（i）不动，沿着行（j）走    for i in range(sy):        for j in range(sx):            result[j:j + p, i:i + q] += mtx[:, col].reshape(block_size, order='F')            weight[j:j + p, i:i + q] += np.ones(block_size)            col += 1    return result / weightdef sparse_encode(image, block_size, codebook, threshold):    blocks = im2col(image, block_size)    return mp.mpalg.bmp(blocks, codebook, threshold)def sparse_decode(coefficients, codebook, image_size, block_size):    blocks = np.dot(codebook, coefficients)    return col2im(blocks, image_size, block_size)if __name__ == '__main__':    image = mh.imread('Lenna.jpg')    image = mh.colors.rgb2gray(image)    image_size = image.shape    block_size = (8, 8)    codebook = dctmtx(block_size[0] * block_size[1])    threshold = 30    coefficients = sparse_encode(image, block_size, codebook, threshold)    reconstructed = sparse_decode(coefficients, codebook, image_size, block_size)    plt.gray()    plt.subplot(121)    plt.title('原始图像')    plt.imshow(image)    plt.subplot(122)    plt.title('稀疏重建')    plt.imshow(reconstructed)    plt.show()

下面是分别设置threshold为10，20和30的运行结果：

稀疏系数设置为10的重建结果

稀疏系数设置为20的重建结果

稀疏系数设置为30的重建结果

可以看到随着稀疏值的增大，重建的的结果会越来越好，但是稀疏度降低。这中间需要一个平衡。

参考资料

1. 匹配追踪算法原理（GitHub）

2. 匹配追踪算法原理（简书）