基于MATLAB：相干成像以及基于相干衍射成像的光学加密仿真
基于MATLAB：相干成像以及基于相干衍射成像的光学加密仿真
本文将介绍基于MATLAB的相干成像以及基于相干衍射成像的光学加密仿真。首先解释相干成像的基本原理，包括振幅扩散函数和振幅传递函数的概念，并通过MATLAB代码展示相干成像的仿真过程。接着扩展到基于相干衍射成像的光学加密应用，介绍最简单的加密方法，并提供相应的MATLAB代码。
由小编的上一篇文章可知，像差对成像系统的脉冲响应有很强的影响。此外，成像系统模型根据物体照明的相干性有不同的脉冲响应。如果照明是空间相干的，脉冲响应被称为振幅扩散函数（或相干扩散函数），系统的频率响应被称为振幅传递函数（或相干传递函数）。如果照明是空间非相干的，脉冲响应被称为点扩散函数，系统的频率响应被称为光学传递函数（OTF），其幅度被称为调制传递函数（MTF）。需要注意的是，波前像差与照明无关，它们仅取决于成像系统的光学元件。然而，它们对图像的影响却依赖于照明的相干性。
相干成像
当光是相干的，成像系统在光学场中是线性的。因此，图像振幅Ui(u,v)是物体振幅Uo(u,v)与振幅扩散函数h(u,v)的卷积，具体表示为：
考虑成像系统具有单位放大率，振幅扩展函数为：
其中 P(x,y)是一般光瞳函数P(x,y)=A(x,y)exp(iφ(x,y))，zi为像距。如果将卷积定理应用于公式1，结果为：
通过这种形式，我们可以清楚地看到，振幅扩散函数的傅里叶谱调制了物体的频谱，从而生成了衍射图像。振幅传递函数H(fx,fy)描述了成像系统如何在频域中传递和调制物体的频谱。
在仿真过程中，物体由三个平行的矩形狭缝组成，如图1(a)所示。我们假定成像系统中的像差是0.05波长的泽尼克离焦（i=4），由此可以计算出广义瞳孔函数。再使用 ft2 函数计算振幅扩散函数，并在图1(b)中显示。我们发现它比物体要小得多，这是线性系统中脉冲响应的典型特征。最后，使用 conv2 函数将物体场与振幅扩散函数进行卷积，形成最终图像场。最终的物体强度如图1(c)所示。
图1 相干成像实例
MATLAB 仿真代码
% 相干成像
clc;close;clear
N = 256; % 网格数
L = 0.1; % 网格长度 [m]
D = 0.07; % 瞳孔直径 [m]
delta = L / N; % 网格间距[m]
wvl = 1e-6; % 波长 [m]
z = 0.25; % 成像距离 [m]
% pupil-plane coordinates
[x, y] = meshgrid((-N/2 : N/2-1)*delta);
[theta, r] = cart2pol(x, y);
%% 计算振幅扩展函数
% 像差
W = 0.05*zernike(4, 2 * r/D, theta);
% 光瞳函数
P = circ(x, y, D) .* exp(1i*2 * pi*W);
% 振幅扩展函数
h = ft2(P, delta);
%% 物体振幅
%物体坐标和h一致（需要傅里叶变换）
delta_u = wvl*z / (N * delta);
[u v] = meshgrid((-N/2 : N/2-1)*delta_u);
% 物体振幅
obj = (rect((u-1.4e-4)/5e-5) + rect(u/5e-5)+ rect((u+1.4e-4)/5e-5)) .* rect(v/2e-4);
%% 相干成像
img = myconv2(obj, h, 1);
figure
subplot(1,3,1)
%物体振幅
imagesc(abs(obj));axis off ;colormap("parula");axis square
subplot(1,3,2)
%振幅扩展函数
imagesc(abs(h));axis off ;colormap("parula");axis square
%成像振幅
subplot(1,3,3)
imagesc(abs(img));axis off ;colormap("parula");axis square
相干衍射成像加密
上文展示了物体振幅通过成像系统相干成像的过程，在此基础上，我们可以扩展到基于相干衍射成像光学加密的应用。光学加密有很多方法，这里我们只展示最简单的方法，就是直接对已知图像做傅里叶变换得到一张加密后的图像，经过成像系统成像后，接收到的图像进行逆傅里叶变换便完成解密得到我们想要的图像。仿真结果展示在图2中。实例中的加密方法，振幅扩散函数等可以自行改变，此外可以将图像联合相位进行加密，在接收端利用GS算法等相位恢复算法对图像进行解密。非相干成像的内容敬请期待下一期。
MATLAB 仿真代码
% 相干成像加密
clc;close;clear
N = 320; % 网格数
L = 0.1; % 网格长度 [m]
D = 0.1; % 瞳孔直径 [m]
delta = L / N; % 网格间距[m]
wvl = 1e-6; % 波长 [m]
z = 0.25; % 成像距离 [m]
% pupil-plane coordinates
[x, y] = meshgrid((-N/2 : N/2-1)*delta);
[theta, r] = cart2pol(x, y);
%% 计算振幅扩展函数
% 像差
W = 0.05*zernike(4, 2 * r/D, theta);
% 光瞳函数
P = circ(x, y, D) .* exp(1i*2 * pi*W);
% 振幅扩展函数
h = ft2(P, delta);
%% 需要传输的图像
%物体坐标和h一致（需要傅里叶变换）
delta_u = wvl*z / (N * delta);
[u v] = meshgrid((-N/2 : N/2-1)*delta_u);
%加密图像
obj=imread("cat.jpg");
obj=double(rgb2gray(obj));%做灰度处理，转为double类型便于计算
objencode=ft2(obj,delta_u);
%% 相干成像并解密
img = myconv2(objencode, h, 1);
imgf=ift2(img,1);
figure
subplot(2,2,1)
%物体振幅
imagesc(abs(obj));axis off ;colormap("gray");axis square
subplot(2,2,2)
%加密后的图像
imagesc(abs(objencode));axis off ;colormap("gray");axis square
%成像振幅
subplot(2,2,3)
imagesc(abs(img));axis off ;colormap("gray");axis square
%解密后的图像
subplot(2,2,4)
imagesc(abs(imgf));axis off ;colormap("gray");axis square