【youcans 的 OpenCV 例程 300 篇】106. 退化图像的逆滤波

6. 退化图像复原

图像复原是对图像退化的过程进行估计，并补偿退化过程造成的失真，以便获得未经退化的原始图像或原始图像的最优估值，从而改善图像质量的一种方法。

典型的图像复原方法是根据图像退化的先验知识建立退化模型，以退化模型为基础采用滤波等手段进行处理，使复原后的图像符合一定的准则，达到改善图像质量的目的。

因此，图像复原是沿着质量降低的逆过程来重现真实的原始图像，通过去模糊函数而去除图像模糊。

6.1 退化图像的逆滤波（Inverse filter）

图像退化表示为退化算子 $H\mathcal{H}$ ，退化函数可以用观察法、试验法或建模法估计，则通过逆滤波就可以直接实现图像复原。用退化图像的傅里叶变换除以退化函数的傅里叶变换，得到原始图像的傅里叶变换估计：
$F^(u,v)=G(u,v)H(u,v)\hat{F}(u,v) = \frac{G(u,v)}{H(u,v)}$

但是，由于实际上退化图像是退化算子与加性噪声项共同作用的结果，由此得到：
$F^(u,v)=F(u,v)+N(u,v)H(u,v)\hat{F}(u,v) = F(u,v) + \frac{N(u,v)}{H(u,v)}$

这表明即使获得退化函数 $H (u, v)$ 的估计，由于噪声项是未知的，因此也不能准确地复原原始图像。

进一步地，如果退化函数为 0 或很小，则噪声项的影响将非常严重（信噪比低）。这时，需要将频率限制到原点附近进行分析，可以减少遇到零值的可能性。

例程 9.20：湍流模糊退化图像的逆滤波

如前所述，通过湍流退化模型可以得到退化图像。使用该退化模型进行逆滤波，退化函数与生成退化图像所用的退化函数相反：
$H(u,v) = e^{-k [(u-M/2)^2+(v-N/2)^2]^{5/6}}$

但是，直接使用退化模型 H(u,v) 逆滤波的结果（D0=full）很差，用理想低通滤波器对退化模型 H(u,v) 在半径 D0 之外截止后，则视觉效果较好。

    # 9.20: 湍流模糊退化图像的逆滤波def turbulenceBlur(img, k=0.001):  # 湍流模糊传递函数: H(u,v) = exp(-k(u^2+v^2)^5/6)M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))radius = (u - M//2)**2 + (v - N//2)**2kernel = np.exp(-k * np.power(radius, 5/6))return kerneldef getDegradedImg(image, Huv, eps):  # 根据退化模型生成退化图像# (1) 傅里叶变换, 中心化fft = np.fft.fft2(image.astype(np.float32))  # 傅里叶变换fftShift = np.fft.fftshift(fft)  # 将低频分量移动到频域图像中心# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数fftShiftFilter = fftShift * Huv  # Guv = Fuv * Huv# (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化invShift = np.fft.ifftshift(fftShiftFilter)  # 将低频分量逆转换回图像四角imgIfft = np.fft.ifft2(invShift)  # 逆傅里叶变换，返回值是复数数组imgDegraded = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return imgDegradeddef ideaLPFilter(img, radius=10):  # 理想低通滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))D = np.sqrt((u - M//2)**2 + (v - N//2)**2)kernel = np.zeros(img.shape[:2], np.float32)kernel[D <= radius] = 1return kerneldef inverseFilter(image, Huv, D0):  # 根据退化模型逆滤波# (1) 傅里叶变换, 中心化fft = np.fft.fft2(image.astype(np.float32))  # 傅里叶变换fftShift = np.fft.fftshift(fft)  # 将低频分量移动到频域图像中心# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数if D0==0:fftShiftFilter = fftShift / Huv  # Guv = Fuv / Huvelse:lpFilter = ideaLPFilter(image, radius=D0)fftShiftFilter = fftShift / Huv * lpFilter  # Guv = Fuv / Huv# (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化invShift = np.fft.ifftshift(fftShiftFilter)  # 将低频分量逆转换回图像四角imgIfft = np.fft.ifft2(invShift)  # 逆傅里叶变换，返回值是复数数组imgRebuild = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return imgRebuild# 读取原始图像img = cv2.imread("../images/Fig0525a.tif", 0)  # flags=0 读取为灰度图像# 生成湍流模糊图像HTurb = turbulenceBlur(img, k=0.0025)imgBlur = np.abs(getDegradedImg(img, HTurb, 0.0))print(imgBlur.max(), imgBlur.min())# # 逆滤波imgRebuild = inverseFilter(imgBlur, HTurb, 480)  # Huv 全滤波器imgRebuild1 = inverseFilter(imgBlur, HTurb, D0=40)  # 在半径 D0 之外 Huv 截止imgRebuild2 = inverseFilter(imgBlur, HTurb, D0=70)imgRebuild3 = inverseFilter(imgBlur, HTurb, D0=100)plt.figure(figsize=(9, 7))plt.subplot(231), plt.title("origin"), plt.axis('off'), plt.imshow(img, 'gray')plt.subplot(232), plt.title("turbulence blur"), plt.axis('off'), plt.imshow(imgBlur, 'gray')plt.subplot(233), plt.title("inverse filter(D0=full)"), plt.axis('off'), plt.imshow(imgRebuild, 'gray')plt.subplot(234), plt.title("inverse filter(D0=40)"), plt.axis('off'), plt.imshow(imgRebuild1, 'gray')plt.subplot(235), plt.title("inverse filter(D0=70)"), plt.axis('off'), plt.imshow(imgRebuild2, 'gray')plt.subplot(236), plt.title("inverse filter(D0=100)"), plt.axis('off'), plt.imshow(imgRebuild3, 'gray')plt.tight_layout()plt.show()