一：图像复原代数方法

图像复原代数方法：根据退化模型，假设具备关于$g$、$H$、$n$的某些先验知识，确定某种最佳准则，寻找原图像$f$的最优估计

（1）无约束最小乘方复原

无约束最小乘方复原：是一种用于恢复受损图像的方法。它基于最小化平方误差的原则，以尽可能接近原始图像为目标。假设我们有一个受损的图像表示为矩阵$Y$，其大小为$m\times n$，其中的元素表示像素灰度值。我们希望通过图像复原算法来估计一个恢复的图像矩阵$X$，大小也为$m\times n$。我们可以将图像复原看作是一个逆问题，其中观测数据是受损的图像矩阵$Y$，而未知量是恢复的图像矩阵$X$。我们假设受损图像$Y$与恢复图像$X$之间存在线性关系，可以用一个称为复原算子$H$的矩阵来表示。那么，无约束最小乘方复原的目标是找到一个估计的恢复图像矩阵$X^{\wedge }$，使得通过复原算子$H$和估计的恢复图像矩阵$X^{\wedge }$计算得到的预测图像矩阵$Y^{\wedge }$与观测图像矩阵Y之间的平方误差最小。数学上可以表示为

$$X^{\wedge }=\mathrm{argmin}\Vert Y-HX{\Vert }^2$$

通过求解上述优化问题，我们可以得到估计的恢复图像矩阵$X^{\wedge }$，它是使得平方误差最小的解。无约束最小乘方复原方法可基于正则化方法、优化算法等进行求解，以获得更好的复原效果

（2）约束复原

约束复原：附加某种约束条件，称为约束复原。有附加条件的极值问题可用拉格朗日乘数法来求解

• 约束条件：$\Vert n{\Vert }^2=\Vert g-H\hat{f}{\Vert }^2$
• 准则函数：$\Vert Q\hat{f}{\Vert }^2$，$Q$为对原图像进行的某一线性运算
• 拉格朗日函数：$J(\hat{f})=\Vert Q\hat{f}{\Vert }^2+\lambda \left(\Vert g-H\hat{f}{\Vert }^2-\Vert n{\Vert }^2\right)$
• 求极小值：$\hat{f}={\left(H^{T}H+\frac{1}{\lambda }{Q}^{T}Q\right)}^{-1}{H}^{T}g$

二：典型图像复原方法

（1）逆滤波复原

A：概述

逆滤波复原：是图像复原中一种常用的方法，旨在通过对图像进行频域上的逆滤波操作，恢复原始图像。假设我们有一个受损的图像表示为矩阵$Y$，大小为$m\times n$，其中的元素表示像素灰度值。我们希望通过逆滤波复原算法来估计一个恢复的图像矩阵$X$，大小也为$m\times n$。逆滤波复原的基本原理是假设图像受损是由于在频域上的卷积模糊操作所引起的。我们将这个模糊操作的效应表示为一个称为模糊核或点扩散函数（Point Spread Function，PSF）的矩阵，记作$H$。那么，在频域上的逆滤波复原可以表示为以下数学表达式

• $F$：离散傅里叶变换
• $F^{-1}$：离散傅里叶逆变换
• $F(Y)$和$F(H)$分别表示对图像矩阵$Y$和模糊核矩阵$H$进行离散傅里叶变换得到的频域表示

$$X^{\wedge }=F^{-1}\frac{F(Y)}{F(H)}$$

在这个逆滤波复原公式中，我们将观测图像矩阵$Y$和模糊核矩阵$H$都转换到频域上，然后通过将它们的频域表示相除来进行逆滤波操作。最后，再将逆滤波后的结果通过离散傅里叶逆变换恢复到空域，得到估计的恢复图像矩阵$X^{\wedge }$

B：程序

如下：对图像进行均值模糊，并进行逆滤波复原

---

matlab实现：

clear,clc,close all;
Image=im2double(rgb2gray(imread('flower.jpg')));
window=15;
[n,m]=size(Image);
n=n+window-1;
m=m+window-1;
h=fspecial('average',window);
BlurI=conv2(h,Image);
BlurandnoiseI=imnoise(BlurI,'salt & pepper',0.001);
figure,imshow(Image),title('Original Image');
figure,imshow(BlurI),title('Blurred Image');
figure,imshow(BlurandnoiseI),title('Blurred Image with noise');
imwrite(BlurI,'Blurred.jpg');
imwrite(BlurandnoiseI,'blurrednoise.jpg');

h1=zeros(n,m);
h1(1:window,1:window)=h;
H=fftshift(fft2(h1));
H(abs(H)<0.0001)=0.01;
M=H.^(-1);
r1=floor(m/2);  r2=floor(n/2);d0=sqrt(m^2+n^2)/20;
for u=1:m
    for v=1:n
        d=sqrt((u-r1)^2+(v-r2)^2);
        if d>d0
            M(v,u)=0;
        end
    end
end

G1=fftshift(fft2(BlurI));
G2=fftshift(fft2(BlurandnoiseI));

f1=ifft2(ifftshift(G1./H));
f2=ifft2(ifftshift(G2./H));
f3=ifft2(ifftshift(G2.*M));

result1=f1(1:n-window+1,1:m-window+1);
result2=f2(1:n-window+1,1:m-window+1);
result3=f3(1:n-window+1,1:m-window+1);

figure,imshow(abs(result1),[]),title('Filtered Image1');
figure,imshow(abs(result2),[]),title('Filtered Image2');
figure,imshow(abs(result3),[]),title('Filtered Image3');

imwrite(abs(result1),'Filteredimage1.jpg');
imwrite(abs(result2),'Filteredimage2.jpg');
imwrite(abs(result3),'Filteredimage3.jpg');


% clear,clc;
% f=[1 1 1 1 1 1 1 1];
% h=[1 1 1];
% f1=[f 0 0];
% h1=[h 0 0 0 0 0 0 0];
% g1=conv(h,f);
% G1=fft(g1);
% H1=fft(h1);
% M=H1.^(-1);
% M(5:10)=0;
% F1=G1./H1;
% f2=ifft(F1);
% F2=G1.*M;
% f3=ifft(F2);

Python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

image = cv2.imread('flower.jpg')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
image = image.astype(np.float64) / 255.0

window = 15
n, m = image.shape
n = n + window - 1
m = m + window - 1

h = np.ones((window, window)) / window ** 2
blur_image = cv2.filter2D(image, -1, h)

blur_noise_image = np.copy(blur_image)
noise = np.random.rand(*blur_noise_image.shape)
salt_pepper = (noise < 0.001)
blur_noise_image[salt_pepper] = 0
blur_noise_image[~salt_pepper] = 1

plt.imshow(image, cmap='gray')
plt.title('Original Image')
plt.show()

plt.imshow(blur_image, cmap='gray')
plt.title('Blurred Image')
plt.show()

plt.imshow(blur_noise_image, cmap='gray')
plt.title('Blurred Image with Noise')
plt.show()

cv2.imwrite('Blurred.jpg', blur_image * 255)
cv2.imwrite('blurrednoise.jpg', blur_noise_image * 255)

h1 = np.zeros((n, m))
h1[:window, :window] = h
H = np.fft.fftshift(np.fft.fft2(h1))
H[np.abs(H) < 0.0001] = 0.01
M = np.power(H, -1)

r1 = m // 2
r2 = n // 2
d0 = np.sqrt(m ** 2 + n ** 2) / 20

for u in range(m):
    for v in range(n):
        d = np.sqrt((u - r1) ** 2 + (v - r2) ** 2)
        if d > d0:
            M[v, u] = 0

G1 = np.fft.fftshift(np.fft.fft2(blur_image))
G2 = np.fft.fftshift(np.fft.fft2(blur_noise_image))

f1 = np.fft.ifft2(np.fft.ifftshift(G1 / H))
f2 = np.fft.ifft2(np.fft.ifftshift(G2 / H))
f3 = np.fft.ifft2(np.fft.ifftshift(G2 * M))

result1 = np.abs(f1[:n-window+1, :m-window+1])
result2 = np.abs(f2[:n-window+1, :m-window+1])
result3 = np.abs(f3[:n-window+1, :m-window+1])

plt.imshow(result1, cmap='gray')
plt.title('Filtered Image 1')
plt.show()

plt.imshow(result2, cmap='gray')
plt.title('Filtered Image 2')
plt.show()

plt.imshow(result3, cmap='gray')
plt.title('Filtered Image 3')
plt.show()

cv2.imwrite('Filteredimage1.jpg', result1 * 255)
cv2.imwrite('Filteredimage2.jpg', result2 * 255)
cv2.imwrite('Filteredimage3.jpg', result3 * 255)

（2）维纳滤波复原

A：概述

维也纳滤波复原：是图像复原中一种常用的方法，旨在通过最小化估计图像与原始图像之间的均方误差，对受损图像进行恢复。假设我们有一个受损的图像表示为矩阵$Y$，大小为$m\times n$，其中的元素表示像素灰度值。我们希望通过维纳滤波复原算法来估计一个恢复的图像矩阵$X$，大小也为$m\times n$。维纳滤波复原的基本原理是将复原问题建模为信号的最小均方误差估计问题，通过估计图像与原始图像之间的信噪比以及受损图像与原始图像之间的频率响应之间的关系来进行复原。数学上，维纳滤波复原可以表示为以下表达式

• $X^{\wedge }$：估计的恢复图像矩阵
• $Y$：受损图像矩阵
• $F$：离散傅里叶变换
• $F^{-1}$：离散傅里叶逆变换
• $H$：系统的频率响应
• $S^2$：原始图像的功率谱密度
• $N$：噪声功率谱密度
  $$X^{\wedge }=F^{-1}[H\ast \frac{F(Y)}{H^2+\frac{S^2}{N}}]$$

在这个表达式中，首先通过对受损图像矩阵Y进行离散傅里叶变换得到频域表示，然后将其与系统频率响应$H$相除。这相当于对受损图像在频域上进行去卷积操作。然后，通过对去卷积结果除以噪声功率谱密度和原始图像的功率谱密度的比例，来校正估计图像的信噪比。最后，将校正后的频域结果通过离散傅里叶逆变换转换回空域，得到估计的恢复图像矩阵$X^{\wedge }$

B：程序

如下：对运动模糊的图像进行维纳滤波

---

matlab实现：

deconvwnr函数是MATLAB图像处理工具箱中用于维纳滤波复原的函数之一。它可以用于对受损图像进行复原，通过最小化估计图像与原始图像之间的均方误差，恢复原始图像。该函数的基本语法如下

J = deconvwnr(I, PSF, NSR)

• I：受损图像。可以是灰度图像或彩色图像。
• PSF：点扩散函数（Point Spread Function），即模糊核。可以是灰度图像、二维数组或函数句柄。PSF描述了图像受到的模糊效应。
• NSR：信噪比（Noise-to-Signal Ratio），即受损图像的信噪比。通常使用估计的噪声方差除以信号方差来计算

clear,clc,close all;
Image=im2double(rgb2gray(imread('flower.jpg')));
subplot(221),imshow(Image),title('Original Image');
  
LEN=21;
THETA=11;
PSF=fspecial('motion', LEN, THETA);
BlurredI=imfilter(Image, PSF, 'conv', 'circular');

noise_mean = 0;
noise_var = 0.0001;
BlurandnoisyI=imnoise(BlurredI, 'gaussian', noise_mean, noise_var);
subplot(222), imshow(BlurandnoisyI),title('Simulate Blur and Noise');
imwrite(BlurandnoisyI,'sbn.jpg');

estimated_nsr = 0;
result1= deconvwnr(BlurandnoisyI, PSF, estimated_nsr);
subplot(223),imshow(result1),title('Restoration Using NSR = 0');
imwrite(result1,'runsr0.jpg');
 
estimated_nsr = noise_var / var(Image(:));
result2 = deconvwnr(BlurandnoisyI, PSF, estimated_nsr);
subplot(224),imshow(result2),title('Restoration Using Estimated NSR');
imwrite(result2,'ruensr.jpg');

Python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取原始图像并转为灰度图像
image = cv2.imread('flower.jpg')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
image = image.astype(np.float64) / 255.0

# 显示原始图像
plt.subplot(221)
plt.imshow(image, cmap='gray')
plt.title('Original Image')

# 创建运动模糊核
LEN = 21
THETA = 11
PSF = cv2.getRotationMatrix2D((LEN/2, LEN/2), THETA, 1)
PSF = np.pad(PSF, ((0, LEN-PSF.shape[0]), (0, LEN-PSF.shape[1])), mode='constant', constant_values=0)

# 对原始图像进行模糊处理
blurred_image = cv2.filter2D(image, -1, PSF, borderType=cv2.BORDER_CONSTANT)

# 添加高斯噪声
noise_mean = 0
noise_var = 0.0001
blurred_noisy_image = blurred_image + np.random.normal(noise_mean, np.sqrt(noise_var), blurred_image.shape)

# 显示模拟的模糊和噪声图像
plt.subplot(222)
plt.imshow(blurred_noisy_image, cmap='gray')
plt.title('Simulated Blur and Noise')
plt.imsave('sbn.jpg', blurred_noisy_image, cmap='gray')

# 估计信噪比为0进行复原
estimated_nsr = 0
result1, _ = cv2.deconvolve(blurred_noisy_image, PSF)
result1 = np.clip(result1, 0, 1)

# 显示使用信噪比0进行复原的结果
plt.subplot(223)
plt.imshow(result1, cmap='gray')
plt.title('Restoration Using NSR = 0')
plt.imsave('runsr0.jpg', result1, cmap='gray')

# 估计信噪比并进行复原
estimated_nsr = noise_var / np.var(image)
result2, _ = cv2.deconvolve(blurred_noisy_image, PSF)
result2 = np.clip(result2, 0, 1)

# 显示使用估计的信噪比进行复原的结果
plt.subplot(224)
plt.imshow(result2, cmap='gray')
plt.title('Restoration Using Estimated NSR')
plt.imsave('ruensr.jpg', result2, cmap='gray')

# 显示图像
plt.show()

（3）等功率谱滤波

A：概述

等功率谱滤波：用于降低图像中的噪声或恢复受损的图像。等功率谱滤波的目标是在尽量保持图像的功率谱形状的前提下，减小噪声或修复图像。假设原始图像为$I(x,y)$，其中$(x,y)$是图像中的像素坐标。我们可以将该图像的二维离散傅里叶变换（DFT）表示为$F(u,v)$，其中$(u,v)$是频域中的频率坐标。等功率谱滤波的思想是，我们可以通过适当地缩放原始图像的频域表示来减小噪声或修复图像。具体而言，我们可以将频域表示$F(u,v)$乘以一个滤波函数$H(u,v)$，以实现等功率谱滤波。滤波函数$H(u,v)$用于控制不同频率成分的增益或衰减。其数学表达式如下，$G(u,v)$是经过等功率谱滤波后的频率表示

$$G(u,v)=F(u,v)\ast H(u,v)$$

为了满足等功率谱的条件，我们希望滤波函数$H(u,v)$满足以下关系

$$|G(u,v){|}^2\approx |F(u,v){|}^2$$

也即经过滤波后的频域表示的功率谱应接近原始频域表示的功率谱。这样可以保持图像的能量分布，减小噪声或修复受损的图像。然后，通过应用逆离散傅里叶变换（IDFT）将滤波后的频域表示$G(u,v)$转换回空域表示，可以得到经过等功率谱滤波的图像

B：程序

如下：对运动模糊加噪声图像进行等功率谱滤波复原

---

matlab实现：

clear,clc,close all;
Image=im2double(rgb2gray(imread('flower.jpg')));
[n,m]=size(Image);
figure,imshow(Image),title('Original Image');

LEN=21;
THETA=11;
PSF=fspecial('motion', LEN, THETA);
BlurredI=conv2(PSF,Image);
figure,imshow(BlurredI),title('motion blur');
imwrite(BlurredI,'motionblur.jpg');

[nh,mh]=size(PSF);
n=n+nh-1;
m=m+mh-1;
noise=imnoise(zeros(n,m),'salt & pepper',0.001);
BlurandnoiseI=BlurredI+noise;
figure,imshow(BlurandnoiseI),title('Blurred Image with noise');  
imwrite(BlurandnoiseI,'motionblurandnoise.jpg');

h1=zeros(n,m);
h1(1:nh,1:mh)=PSF;
H=fftshift(fft2(h1));
K=sum(noise(:).^2)/sum(Image(:).^2);
M=(1./(abs(H).^2+K)).^0.5;

G=fftshift(fft2(BlurandnoiseI));
f=ifft2(ifftshift(G.*M));
result=f(1:n-nh+1,1:m-mh+1);
figure,imshow(abs(result)),title('Filtered Image');
imwrite(abs(result),'denggonglvpuimage.jpg');

Python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像并转为灰度图
image = cv2.imread('flower.jpg')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
image = image.astype(np.float64) / 255.0

# 显示原始图像
plt.imshow(image, cmap='gray')
plt.title('Original Image')
plt.show()

LEN = 21
THETA = 11

# 创建运动模糊点扩散函数（PSF）
PSF = cv2.getMotionKernel(LEN, THETA)
PSF = PSF / np.sum(PSF)  # 归一化

# 进行运动模糊处理
blurred_image = cv2.filter2D(image, -1, PSF)

# 显示模糊图像
plt.imshow(blurred_image, cmap='gray')
plt.title('Motion Blur')
plt.show()

# 保存模糊图像
cv2.imwrite('motionblur.jpg', blurred_image * 255.0)

nh, mh = PSF.shape
n, m = image.shape
n = n + nh - 1
m = m + mh - 1

# 添加椒盐噪声
noise = np.random.rand(n, m)
noise = np.where(noise < 0.001, 1, 0)
blurred_noisy_image = blurred_image + noise

# 显示带噪声的模糊图像
plt.imshow(blurred_noisy_image, cmap='gray')
plt.title('Blurred Image with Noise')
plt.show()

# 保存带噪声的模糊图像
cv2.imwrite('motionblurandnoise.jpg', blurred_noisy_image * 255.0)

h1 = np.zeros((n, m))
h1[0:nh, 0:mh] = PSF

# 计算滤波器的频域表示
H = np.fft.fftshift(np.fft.fft2(h1))

# 计算噪声方差和信号方差的比例
K = np.sum(noise**2) / np.sum(image**2)

# 计算等功率谱滤波器
M = (1.0 / (np.abs(H)**2 + K))**0.5

# 对带噪声的模糊图像进行频域滤波
G = np.fft.fftshift(np.fft.fft2(blurred_noisy_image))
f = np.fft.ifft2(np.fft.ifftshift(G * M))

# 提取滤波结果的有效区域
result = np.abs(f[0:n-nh+1, 0:m-mh+1])

# 显示滤波后的图像
plt.imshow(result, cmap='gray')
plt.title('Filtered Image')
plt.show()

# 保存滤波后的图像
cv2.imwrite('denggonglvpuimage.jpg', result * 255.0)

（4）几何均值滤波

集合均值滤波：在图像复原中，几何均值滤波是一种常用的滤波方法，用于降低图像中的噪声或改善图像质量。几何均值滤波的主要思想是通过计算像素周围邻域的几何平均值来平滑图像，从而抑制噪声并保留图像的边缘信息。假设原始图像为$I(x,y)$，其中($(x,y)$是图像中的像素坐标。几何均值滤波器的邻域大小为$n\times n$，以中心像素为基准，我们可以定义一个邻域矩阵为$N(x,y$。几何均值滤波的计算步骤如下

• 对于图像中的每个像素点$(x,y)$，将邻域矩阵$N(x,y)$中的所有像素值相乘得到乘积$P$。也即$P=\prod N(i,j)$
• 计算领域矩阵$N(x,y)$中像素值的几何平均值$A$。$A=P^{\frac{1}{n^2}}$
• 将几何平均值$A$作为滤波后的像素值，用于替换原始图像中的像素值$I(x,y)$

几何均值滤波通过计算像素周围邻域的几何平均值来实现平滑处理。由于几何平均值对噪声像素具有较强的抑制作用，因此几何均值滤波常用于去除椒盐噪声等图像中的随机噪声。它能够保持图像的整体亮度和结构特征，但可能导致图像细节的模糊。因此，在应用几何均值滤波时需要权衡噪声抑制和细节保留之间的平衡

（5）约束最小二乘法滤波

A：概述

约束最小二乘法滤波：图像复原中的约束最小二乘方滤波是一种常用的滤波方法，用于恢复受损的图像。该方法基于最小二乘方准则，通过约束滤波器的响应以尽量满足先验信息和图像特性，从而减小噪声和估计图像的失真。假设原始图像为$I(x,y)$，其中$(x,y)$是图像中的像素坐标。我们的目标是通过滤波操作来估计原始图像。首先，我们引入一个约束函数$C$，用于限制滤波器的响应。约束函数$C$可以是对滤波器的频率响应进行限制的函数，例如限制滤波器的幅度响应或相位响应。滤波器的频率响应为$H(u,v)$，其中$(u,v)$是频率域中的频率坐标。通过应用约束函数$C$，我们可以得到约束后的滤波器响应

$$H_c(u,v)=C[H(u,v)]$$
然后，我们通过最小化以下约束最小二乘方问题来估计滤波后的图像

$$\arg \underset{F}{\min }\sum \sum {\left|G(u,v)-F(u,v)\cdot H_c(u,v)\right|}^2$$

其中，$G(u,v)$是受损图像的频率表示，$F(u,v)$是滤波后的图像的频率表示。约束最小二乘方滤波的目标是找到滤波后的图像$F(u,v)$，使得滤波后的图像与受损图像之间的差异最小。最小二乘方问题的解可以通过将其转化为频率域中的点对点相除来获得

$$F(u,v)=\frac{G(u,v)}{H_c(u,v)}$$

通过将滤波后的图像的频率表示$F(u,v)$转换回空域表示，可以得到估计的复原图像

B：程序

如下：对模糊的图像进行约束最小二乘方滤波

---

matlab实现：

deconvreg函数用于执行正则化盲去卷积（regularized blind deconvolution）操作。这个函数可以用于估计图像的退化过程（如模糊或噪声）并恢复原始图像。其语法格式如下

[J, LAMBDA] = deconvreg(I, PSF, K)
[J, LAMBDA] = deconvreg(I, PSF, K, OPTIONS)

输入参数

• I：观测图像，即受退化影响的图像。
• PSF：点扩散函数（Point Spread Function），表示退化过程中的模糊过程。
• K：正则化参数，用于平衡图像恢复和噪声增强之间的权衡

可选参数

• OPTIONS：一个结构体，包含额外的选项设置，如迭代次数、停止准则等

输出参数

• J：恢复的图像，即通过盲去卷积估计得到的原始图像。
• LAMBDA：正则化参数的最佳值，用于平衡去卷积和噪声增强

clear,clc,close all;
Image=im2double(rgb2gray(imread('flower.jpg')));
window=15;
[N,M]=size(Image);
N=N+window-1;
M=M+window-1;
h=fspecial('average',window);
BlurI=conv2(h,Image);
sigma=0.001;
miun=0;
nn=M*N*(sigma+miun*miun);
BlurandnoiseI=imnoise(BlurI,'gaussian',miun,sigma);
figure,imshow(BlurandnoiseI),title('Blurred Image with noise');
imwrite(BlurandnoiseI,'Blurrednoiseimage.jpg');

h1=zeros(N,M);
h1(1:window,1:window)=h;
H=fftshift(fft2(h1));

lap=[0 1 0;1 -4 1;0 1 0];
L=zeros(N,M);
L(1:3,1:3)=lap;
L=fftshift(fft2(L));

G=fftshift(fft2(BlurandnoiseI));
gama=0.3;
step=0.01;
alpha=nn*0.001;
flag=true;
while flag 
    MH=conj(H)./(abs(H).^2+gama*(abs(L).^2));
    F=G.*MH;
    E=G-H.*F;
    E=abs(ifft2(ifftshift(E)));
    ee=sum(E(:).^2);
    if ee<nn-alpha
        gama=gama+step;
    elseif ee>nn+alpha
        gama=gama-step;
    else
        flag=false;
    end
end
MH=conj(H)./(abs(H).^2+gama*(abs(L).^2));
f=ifft2(ifftshift(G.*MH));
result=f(1:N-window+1,1:M-window+1);
[J, LAGRA]=deconvreg(BlurandnoiseI,h,nn);
figure,imshow(J,[]);

figure,imshow(abs(result),[]),title('Filtered Image');
imwrite(abs(result),'LSFilteredimage.jpg');

python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

image = cv2.imread('flower.jpg')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
image = image.astype(np.float64) / 255.0

window = 15
N, M = image.shape
N = N + window - 1
M = M + window - 1

h = cv2.getGaussianKernel(window, 0)
h1 = np.zeros((N, M))
h1[:window, :window] = h

BlurI = cv2.filter2D(image, -1, h)
sigma = 0.001
miun = 0
nn = M * N * (sigma + miun * miun)

BlurandnoiseI = cv2.randn(BlurI, miun, sigma)
BlurandnoiseI = BlurandnoiseI.astype(np.float64) / 255.0

plt.imshow(BlurandnoiseI, cmap='gray')
plt.title('Blurred Image with noise')
plt.show()

H = np.fft.fftshift(np.fft.fft2(h1))

lap = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])
L = np.zeros((N, M))
L[:3, :3] = lap
L = np.fft.fftshift(np.fft.fft2(L))

G = np.fft.fftshift(np.fft.fft2(BlurandnoiseI))
gama = 0.3
step = 0.01
alpha = nn * 0.001
flag = True

while flag:
    MH = np.conj(H) / (np.abs(H) ** 2 + gama * (np.abs(L) ** 2))
    F = G * MH
    E = G - H * F
    E = np.abs(np.fft.ifft2(np.fft.ifftshift(E)))
    ee = np.sum(E ** 2)

    if ee < nn - alpha:
        gama = gama + step
    elif ee > nn + alpha:
        gama = gama - step
    else:
        flag = False

MH = np.conj(H) / (np.abs(H) ** 2 + gama * (np.abs(L) ** 2))
f = np.fft.ifft2(np.fft.ifftshift(G * MH))
result = f[:N-window+1, :M-window+1]

plt.imshow(np.abs(result), cmap='gray')
plt.title('Filtered Image')
plt.show()

（6）Richardson-Lucy算法

A：概述

Richardson-Lucy算法：RL算法是一种迭代盲去卷积算法，用于图像复原。它被广泛用于恢复被模糊和受损的图像。假设我们有一个受模糊和受损的图像$I$，我们的目标是估计原始图像$J$，其中$J$是我们要恢复的图像。首先，我们定义一个点扩散函数（Point Spread Function）H，表示模糊过程中的系统响应。RL算法迭代步骤如下

• 初始化：将估计的图像J初始化为与输入图像I相同的值或使用其他合适的初始化方法。
• 步骤一：计算反卷积因子D。这可以通过将输入图像I与估计的图像J进行卷积来得到，即 D = I ⊗ J。
• 步骤二：计算更新系数E。这可以通过将反卷积因子D与模糊函数H进行卷积来得到，即 E = D ⊗ H。
• 步骤三：计算更新的图像估计J’。这可以通过将估计的图像J与更新系数E逐元素相除来得到，即 J’ = J ⊘ E
• 重复步骤二和步骤三，直到满足停止准则（例如达到最大迭代次数或估计的图像收敛）

B：程序

如下：对模糊的图像进行RL算法复原

---

matlab实现：

deconvlucy函数是用于执行Lucy-Richardson算法的函数。Lucy-Richardson算法是一种迭代盲去卷积算法，用于图像复原。其语法格式如下

[J, PSF] = deconvlucy(I, PSF, NUMIT)
[J, PSF] = deconvlucy(I, PSF, NUMIT, DAMPAR)
[J, PSF] = deconvlucy(I, PSF, NUMIT, DAMPAR, INITPSF)
[J, PSF] = deconvlucy(I, PSF, NUMIT, DAMPAR, INITPSF, OPTIONS)

输入参数：

• I：观测图像，即受退化影响的图像。
• PSF：点扩散函数（Point Spread Function），表示退化过程中的模糊过程。
• NUMIT：迭代次数，指定Lucy-Richardson算法的迭代次数。

可选参数：

• DAMPAR：阻尼参数（Damping Parameter），用于控制算法的收敛性和稳定性。
• INITPSF：初始化的点扩散函数，用于改进算法的初始估计。
• OPTIONS：一个结构体，包含其他选项设置，如停止准则、正则化等。

输出参数：

• J：恢复的图像，即通过Lucy-Richardson算法估计得到的原始图像。
• PSF：估计的点扩散函数

clear,clc,close all;
I=im2double(rgb2gray(imread('flower.jpg')));
figure,imshow(I),title('原图像');
PSF=fspecial('gaussian',7,10);%产生一个高斯低通滤波器，模板尺寸为[7 7],滤波器的标准差为10
V=0.0001;%高斯加性噪声的标准差
IF1=imfilter(I,PSF);%原图像通过高斯低通滤波器
BlurredNoisy=imnoise(IF1,'gaussian',0,V);%加入高斯噪声
figure,imshow(BlurredNoisy),title('高斯模糊加噪声图像');

WT=zeros(size(I));%产生权重矩阵
WT(5:end-1,5:end-4)=1;

%使用不同的参数进行复原
J1=deconvlucy(BlurredNoisy,PSF);
J2=deconvlucy(BlurredNoisy,PSF,50,sqrt(V));
J3=deconvlucy(BlurredNoisy,PSF,100,sqrt(V),WT);

figure,imshow(J1),title('10次迭代');
figure,imshow(J2),title('50次迭代');
figure,imshow(J3),title('100次迭代');

python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像
I = cv2.imread('flower.jpg')
I = cv2.cvtColor(I, cv2.COLOR_BGR2GRAY)
I = I.astype(np.float64) / 255.0

# 显示原图像
plt.imshow(I, cmap='gray')
plt.title('原图像')
plt.show()

# 创建高斯低通滤波器
PSF = cv2.getGaussianKernel(7, 10)
PSF = np.outer(PSF, PSF)

# 添加高斯噪声
V = 0.0001
IF1 = cv2.filter2D(I, -1, PSF)
BlurredNoisy = cv2.randn(IF1, 0, V)
BlurredNoisy = cv2.add(IF1, BlurredNoisy)

# 显示高斯模糊加噪声图像
plt.imshow(BlurredNoisy, cmap='gray')
plt.title('高斯模糊加噪声图像')
plt.show()

# 创建权重矩阵
WT = np.zeros_like(I)
WT[4:, 4:-3] = 1

# 使用不同参数进行复原
J1 = cv2.deconvolve(BlurredNoisy, PSF)[0]
J2 = cv2.deconvolve(BlurredNoisy, PSF, iterations=50)[0]
J3 = cv2.deconvolve(BlurredNoisy, PSF, iterations=100, reg=np.sqrt(V), weight=WT)[0]

# 显示复原结果
plt.imshow(J1, cmap='gray')
plt.title('10次迭代')
plt.show()

plt.imshow(J2, cmap='gray')
plt.title('50次迭代')
plt.show()

plt.imshow(J3, cmap='gray')
plt.title('100次迭代')
plt.show()