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

空间域上，可以和频域一样进行卷积逆滤波操作。其方法是展开图像为一列，构建卷积模板矩阵，这样卷积操作就变成了矩阵乘法。我们可以用最小二乘法来，已知卷积图像和卷积模板来求出原始图像。空间域最小二乘逆滤波是病态问题，缺点是卷积矩阵非常稀疏和巨大，模非常小，一般需要进行约束。

我们现在有一个均值卷积模板blurKernel，对图像X进行卷积，生成一个模糊图像B，可以写成blurKernel∗X+N=B，其中*是卷积操作,N是加性噪声。展开后我们得到卷积矩阵A，以及B图像的展开Bravel，那么有AXravel+Nravel=Bravel→X′ravel=(ATA+λI)−1ATBravel

假设B是4x4矩阵，比如⎡⎣⎢⎢⎢15913261014371115481216⎤⎦⎥⎥⎥展开后为Bravel=[1,2,3,4,5…16]T；假设A的卷积是3x3的均值模板，不考虑边界，A是6x16矩阵，那么A=⎡⎣⎢1/9......1/91/901/91/91/901/91/91/90000⎤⎦⎥

可以想象到，这些矩阵非常稀疏和巨大。我们可以用稀疏矩阵库来进行计算，其中进行矩阵求逆是最为耗时耗力的一步。

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)逆滤波图像