亚像素级错位图像修复——利用傅里叶变换时移特性

EthanEntropy

已于 2023-04-11 01:27:35 修改

阅读量952

点赞数 2

分类专栏：信号与系统文章标签： matlab

于 2023-04-11 01:23:46 首次发布

本文链接：https://blog.csdn.net/qq_64357371/article/details/130073389

版权

信号与系统专栏收录该内容

1 篇文章 0 订阅

订阅专栏

前言

文章内容：通过matlab实现读取图片并添加亚像素级行间错位，并使用信号与系统知识恢复图像，使用PSNR（峰值信噪比）指标评估恢复效果。

一、图像错位

图像隔行亚像素级错位意即图像矩阵每行间都具有非整数像素的错位量，这种错位无法直接对图像矩阵操作实现，需要将图像每行的值类比时域的离散时间信号值，对其进行傅里叶变换后通过乘上线性相位的方式，让其发生“时移”（即错位），通过调整线性相位的大小可以实现非整数错位，即亚像素级错位。

具体实现思路：读取bmp类型的Lena图片，将其以double类型存储，存为512*512的方阵。对图像每行进行傅里叶变换后，依次在频域添加随机线性相位，随机量介于0-30个像素之间，再对其进行傅里叶逆变换得到错位的图像。

matlab语言编写的代码如下：

%图像亚像素平移在空间域无法实现（因为空间域里图像平移的最小单位就是1个像素）但是可以在频率域里通过加线性相位来实现
g = im2double(imread('lena.bmp'));%将图像转为double类型
f=g;
for i=1:1:512
r(i,:)=fftshift(fft(f(i,:)));
end
for i=2:1:512
deltat=(2*rand()-1)*30;
for w=1:1:512
r(i,w)=r(i,w)*exp(-1j*(deltat*(w-256)*2*pi/512));
end
end
for i=1:1:512
g(i,:)=ifft(r(i,:));
end

效果如下图：

二、图像复原

相邻两行图像空间域形状相近，其轮廓由傅里叶变换后的低频分量决定，通过使用atan2()函数求各行相位并让相邻两行作差，可以求出每行相对上一行的错位量。

为了还原图像，只需要考虑低频分量（该项目选用基波分量）对应的线性相位计算得到错位量，在频域通过减去该线性相位的方式复原图像。具体实现方法如下：

利用atan2()函数求相邻两行的基波分量的相位，作差得到线性相位因子delta_w，反映第二行错位量。通过在频域让某行减去线性相位因子，可以将图像还原，再对每行逆变换即得到复原后的图像。代码如下：

%图像恢复
repic(1,:)=f(1,:);
for i=1:1:511
fai2=atan2(imag(r(i+1,257)),real(r(i+1,257)));
fai1=atan2(imag(r(i,257)),real(r(i,257)));
delta_w = fai2-fai1;
time_shift(i,1)=delta_w*(256/pi)/(-1);
for k=1:512
r(i+1,k)=r(i+1,k).*exp(-1i*delta_w*(k-256));
end
repic(i+1,:)=ifft(r(i+1,:));
end
repic = abs(repic);

每行错位量可通过公式计算，将其可视化：

数值大小反映错位像素点数多少。

复原后图像和错位图对比：

可以看出效果非常好。

三、评估效果

PSNR是“Peak Signal to Noise Ratio”的缩写，即峰值信噪比。其值反映两张图像结构相似性，是图像压缩和复原领域的典型指标。峰值信噪比PSNR介于20dB到30dB之间时，人眼可以察觉出图像差异，大于30dB而小于50dB时，人眼难以察觉图像差异，接近50dB（在50dB附近）时，代表两图误差很小，若PSNR更大，则表明两图几乎无任何误差，完全可以忽略。其计算公式如下：

公式中的MSE表示两个图像中每个相同位置像素值相减平方再求和然后求平均值，即经过平均化操作的err值。表的的是两幅图每个位置上像素值的平均差异，MSE越小，表示相似度越高，对应的PSNR越大。由于MSE对平移敏感的特性，PSNR也对错位图复原效果较敏感，是很好的质量评估指标。代码如下：

%PSNR 峰值信噪比
if sum(sum(img1-img2)) == 0
error('Those pictures are the same');
end
MAX=1; %图像有多少灰度级
if (max(max(img1))-min(min(img1))) ~= 0
img1 = (img1-min(min(img1)))./(max(max(img1))-min(min(img1)));
end
if (max(max(img1))-min(min(img1))) ~= 0
img2 = (img2-min(min(img2)))./(max(max(img2))-min(min(img2)));
end 
MSE=sum(sum((img1-img2).^2))/(512*512); %图片像素
PSNR=20*log10(MAX/sqrt(MSE)); %峰值信噪比

经计算，该次测试中PSNR为301.2313147989404，非常大，表明图像复原效果非常好。

总结

本文探讨了用信号与系统的方法复原具有隔行亚像素级错位的灰度图的具体案例和实现流程。这种方法简单可靠、效果良好。

在实际工程领域，隔行亚像素错位图像的复原也有现实案例：扫描雷达成像时由于雷达设备的机械振动和摩擦等设备因素，可能会出现成像图出现行间错位的情况。在解决这类问题的时候，用这种方法进行后期的图像复原处理从时间和成本上来看是很不错的方法。不过，由于在现实场景中存在各种干扰因素，其复原质量不会像本文案例这么理想，读者可思考更多效果更好、鲁棒性更强的方法。

附录

完整程序代码如下：

clear all; 
%图像亚像素平移在空间域无法实现（因为空间域里图像平移的最小单位就是1个像素）但是可以在频率域里通过加线性相位来实现
g = im2double(imread('lena.bmp'));%将图像转为double类型
f=g;
for i=1:1:512
r(i,:)=fftshift(fft(f(i,:)));
end
for i=2:1:512
    deltat=(2*rand()-1)*30;
    for w=1:1:512
    r(i,w)=r(i,w)*exp(-1j*(deltat*(w-256)*2*pi/512));
    end
end
for i=1:1:512
g(i,:)=ifft(r(i,:));
end

%图像恢复
repic(1,:)=f(1,:);
for i=1:1:511
    fai2=atan2(imag(r(i+1,257)),real(r(i+1,257)));
    fai1=atan2(imag(r(i,257)),real(r(i,257)));
    delta_w = fai2-fai1;
    time_shift(i,1)=delta_w*(256/pi)/(-1);
        for k=1:512
            r(i+1,k)=r(i+1,k).*exp(-1i*delta_w*(k-256));
        end
    
    repic(i+1,:)=ifft(r(i+1,:));
end
repic = abs(repic);

figure(1);
subplot(1,2,1);
imshow(abs(f));
title('Reference image, f(x,y)')
subplot(1,2,2);
imshow(abs(g));
title('Shifted image, g(x,y)') 
figure(2);
imshow(abs(repic));
title('Fixed image, repic(x,y)') 


figure(3)
time_shift=time_shift';
plot(time_shift(1,:));



%图像恢复质量评估
img1 = f;
img2 = repic;

%误差平方和
err=0;
for i=1:512
    for j=1:512
        err=err+(img1(i,j)-img2(i,j))^2;
    end
end

%PSNR 峰值信噪比
if sum(sum(img1-img2)) == 0
    error('Those pictures are the same');
end
MAX=1;          %图像有多少灰度级
% 归一化
if (max(max(img1))-min(min(img1))) ~= 0
    img1 = (img1-min(min(img1)))./(max(max(img1))-min(min(img1)));
end
if (max(max(img1))-min(min(img1))) ~= 0
    img2 = (img2-min(min(img2)))./(max(max(img2))-min(min(img2)));
end 
MSE=sum(sum((img1-img2).^2))/(512*512);     %图片像素
PSNR=20*log10(MAX/sqrt(MSE));           %峰值信噪比