【图像处理】时域最小二乘逆滤波

空间域上，可以和频域一样进行卷积逆滤波操作。其方法是展开图像为一列，构建卷积模板矩阵，这样卷积操作就变成了矩阵乘法。我们可以用最小二乘法来，已知卷积图像和卷积模板来求出原始图像。空间域最小二乘逆滤波是病态问题，缺点是卷积矩阵非常稀疏和巨大，模非常小，一般需要进行约束。
我们现在有一个均值卷积模板blurKernel，对图像X进行卷积，生成一个模糊图像B，可以写成$blurKernel*X+N=B$，其中*是卷积操作,N是加性噪声。展开后我们得到卷积矩阵A，以及B图像的展开$B_{ravel}$，那么有

$$AX_{ravel}+N_{ravel}=B_{ravel} \\ \rightarrow X_{ravel}' = (A^TA+\lambda I)^{-1}A^TB_{ravel}$$
假设B是4x4矩阵，比如 $\left[ \begin{matrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 \end{matrix} \right]$展开后为 $B_{ravel}=[1, 2, 3, 4, 5 … 16]^T$；假设A的卷积是3x3的均值模板，不考虑边界，A是6x16矩阵，那么 $A=\left[ \begin{matrix} 1/9 & 1/9 & 1/9 & 0 & 1/9 & 1/9 & 1/9 &0 & 1/9 & 1/9 & 1/9 &0 &0 &0 &0 \\ ... \\ ... \end{matrix} \right]$
可以想象到，这些矩阵非常稀疏和巨大。我们可以用稀疏矩阵库来进行计算，其中进行矩阵求逆是最为耗时耗力的一步。

import numpy as np
import os, string
from matplotlib import pyplot as plt
import scipy as sp
import cv2
img = cv2.imread('camera.jpg')
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
img = cv2.resize(img,(0,0),fx=0.5,fy=0.5)
print img.shape
from __future__ import division
blurKernel = np.ones((3,3))/9

from scipy import signal#warning
blurImg = signal.convolve2d(img, blurKernel, 'same','symm')

import itertools
from scipy.sparse import csc_matrix,lil_matrix
from scipy.sparse import linalg as sppl
import numpy.linalg
A=lil_matrix(((img.shape[0]-3)*(img.shape[1]-3),img.shape[0]*img.shape[1]))
print A.shape
ind = -1
for j in xrange(0,np.int32(img.shape[1])-3):
    for i in xrange(0,np.int32(img.shape[0])-3):
        ind +=1
        #index=[]
        ravel = []
        for e in itertools.product(range(j,j+3),range(i,i+3)):
            #index.append(e)
            ravel.append(e[1]+e[0]*img.shape[0])
        A[ind,[ravel]]=1./9.
 A=csc_matrix(A) 

from scipy import sparse
reg = 0.001
X = A.T.dot(A)
X=X+reg*sparse.eye(X.shape[0])
X=csc_matrix.dot(csc_matrix(Xinv), A.T)
padBlur = csc_matrix(blurImg[0:np.int32(img.shape[0])-3,0:np.int32(img.shape[1])-3].ravel()).T
print X.shape,padBlur.shape
X=csc_matrix.dot(X,padBlur)
res=X.todense()
res = res.reshape(blurImg.shape)
fig = plt.figure()
ax = fig.add_subplot(131)  
ax.imshow(img,cmap='gray')
ax = fig.add_subplot(132)  
ax.imshow(blurImg,cmap='gray')
ax = fig.add_subplot(133)  
ax.imshow(res,cmap='gray')
plt.show()

计算花费半个小时。相比频域逆滤波，时间和存储空间代价都很大，一般不推荐直接的空间域逆滤波。可以采用迭代优化的方法。

Fig. (a)原图 (b)均值滤波图像 (c)逆滤波图像