二维图像的离轴数字全息及重建
本文主要参考<衍射计算及数字全息>
[
1
−
3
]
^{[1-3]}
[1−3]一书中相关内容,若有不足的地方,恳请批评指正。
全息摄像主要经历光学全息,数字全息,和计算全息这样几个过程,其基本原理都是通过干涉纪录及衍射重建。当两束光波进行干涉时,纪录物体的振幅和相位信息,通过衍射计算,重构物体的振幅和相位信息。原则上,全息面任一一个点,或者一小块都包含原物体的所有信息,因而,当全息图损坏时,可以通过一小块物体重构原物体。若已知在物平面
y
0
y_{0}
y0,实际模拟运算时,可另物平面为一幅图像,物平面经过衍射距离
d
d
d后,在衍射观测面
y
1
y_{1}
y1的复振幅分布可以以下公式获得。
U
(
p
Δ
x
,
q
Δ
y
)
=
exp
(
j
k
d
)
j
λ
d
exp
[
j
k
2
d
(
p
Δ
x
2
)
+
q
Δ
y
2
]
×
ϝ
[
U
0
(
m
Δ
x
0
,
n
Δ
y
0
)
×
exp
(
j
k
2
d
(
m
Δ
x
0
)
2
,
(
n
Δ
y
0
)
2
)
]
p
Δ
x
λ
d
,
p
Δ
y
λ
d
U(p\Delta x,q\Delta y) = \frac{\exp(jkd)}{j\lambda d}\exp[\frac{jk}{2d}(p\Delta x^2)+q\Delta y^2]\\\times \digamma[U_0(m\Delta x_0,n\Delta y_0)\times\exp(\frac{jk}{2d}(m\Delta x_0)^2,(n\Delta y_0)^2)]_{\frac{p\Delta x}{\lambda d},\frac{p\Delta y}{\lambda d}}
U(pΔx,qΔy)=jλdexp(jkd)exp[2djk(pΔx2)+qΔy2]×ϝ[U0(mΔx0,nΔy0)×exp(2djk(mΔx0)2,(nΔy0)2)]λdpΔx,λdpΔy
具体实现公式1的过程如下:
(1)物平面函数
U
0
(
m
Δ
x
0
,
n
Δ
y
0
)
与
二
次
相
位
因
子
exp
(
j
k
2
d
(
m
Δ
x
0
)
2
,
(
n
Δ
y
0
)
2
)
]
U_0(m\Delta x_0,n\Delta y_0)与二次相位因子\exp(\frac{jk}{2d}(m\Delta x_0)^2,(n\Delta y_0)^2)]
U0(mΔx0,nΔy0)与二次相位因子exp(2djk(mΔx0)2,(nΔy0)2)]相乘,对相乘的结果做傅里叶变换;
(2)对上述的结果与
U
(
p
Δ
x
,
q
Δ
y
)
=
exp
(
j
k
d
)
j
λ
d
exp
[
j
k
2
d
(
p
Δ
x
2
)
+
q
Δ
y
2
]
U(p\Delta x,q\Delta y) = \frac{\exp(jkd)}{j\lambda d}\exp[\frac{jk}{2d}(p\Delta x^2)+q\Delta y^2]
U(pΔx,qΔy)=jλdexp(jkd)exp[2djk(pΔx2)+qΔy2]相乘,即可得到衍射观测面的光波厂复振幅。因为整个计算过程只进行了一次快速傅里叶变换,因而被称为一次傅里叶变换(Singe-FFT)
物面观测面衍射尺寸分析:
之所以要分析物面尺寸关系,是因为不论是数字全息还是计算全息,纪录的相机如CCD,重建的显示设备如空间光调制器LCOS,其靶面尺寸都是有限的。因CCD或者LCOS的靶面尺寸是固定的,因而通过该尺寸反推物面尺寸。选用的CCD的像素尺寸为0.00465mm,像素数为1024则,靶面尺寸为
L
L
L =0.00465×1024;采样数为
N
N
N , 根据S-FFT物平面
L
0
L_0
L0和衍射观测面(观察面)的关系,可以
L
0
=
N
λ
d
L
L_0=\frac{N\lambda d}{L}
L0=LNλd
在衍射重构时,全息面作为物平面; MATLAB模拟程序如下
% h ——波长(mm); I ——数字全息图;
% L ——全息图宽度(mm); z0——记录全息图的距离(mm);
XRGB=imread('baboon.jpg');
X0=rgb2gray(XRGB);
figure,imshow(X0,[]);
[M0,N0]=size(X0);
N1=min(M0,N0);
N=1024; %模拟形成的全息图取样数
X1=imresize(X0,N/4/N1);
[M1,N1]=size(X1);
X=zeros(N,N);
X(N/2-M1/2+1:N/2+M1/2,N/2-N1/2+1:N/2+N1/2)=X1(1:M1,1:N1);
h=0.632e-3; %波长(mm), 可按需要修改
k=2*pi/h;
pix=0.00465 ;%CCD像素宽度(mm), 可按需要修改
L=N*pix; %CCD宽度(mm)
z0=1000; %衍射距离(mm), 可按需要修改
L0=h*N*z0/L; %物平面宽度(mm)
Y=double(X);
a=ones(N,N);
b=rand(N,N)*2*pi;
f=Y.*exp(i.*b); %叠加随机相位噪声,形成振幅正比于图像的初始场复振幅
figstr=strcat('初始物平面宽度=',num2str(L0),'mm');
figure,imshow(X,[]),colormap(gray); xlabel(figstr);title('物平面图像');
%---------------菲涅耳衍射的S-FFT计算开始
n=1:N;
x=-L0/2+L0/N*(n-1);
y=x;
[yy,xx] = meshgrid(y,x); %网格化数据平面
Fresnel=exp(i*k/2/z0*(xx.^2+yy.^2));
f2=f.*Fresnel;
Uf=fft2(f2,N,N);
Uf=fftshift(Uf);
x=-L/2+L/N*(n-1);%CCD宽度取样(mm)
y=x;
[yy,xx] = meshgrid(y,x);
phase=exp(i*k*z0)/(i*h*z0)*exp(i*k/2/z0*(xx.^2+yy.^2));%菲涅耳衍射积分前方的相位因子
Uf=Uf.*phase;
%---------------S-FFT计算结束
figstr=strcat('模拟CCD宽度=',num2str(L),'mm');
figure,imshow(abs(Uf),[]),colormap(gray); xlabel(figstr);title('到达CCD平面的物光振幅分布');
%---------------形成0-255灰度级的数字全息图
fex=N/L;
Qx=(4-2.5)*L0/8/z0; %按照优化设计定义参考光方向余弦
Qy=Qx
x=[-L/2:L/N:L/2-L/N];
y=x;
[X,Y]=meshgrid(x,y);
Ar=max(max(abs(Uf))); %按物光场振幅最大值定义参考光振幅
Ur=Ar*exp(i*k*(X.*Qx+Y.*Qy));%参考光复振幅
Uh=Ur+Uf; %物光与参考光干涉
Wh=Uh.*conj(Uh); %干涉场强度
Imax=max(max(Wh));
I=uint8(Wh./Imax*255); %形成0-255灰度级的数字全息图
imwrite(I,'D:\Ih.tif'); %形成数字全息图文件
figstr=strcat('全息图宽度=',num2str(L),'mm');
figure,imshow(Ih,[]),colormap(gray);xlabel(figstr);title('模拟形成的数字全息图');
%基于S-FFT的全息重建
f0 = double(Ih);
[N1,N2]=size(f0);
N=min(N1,N2);
h=0.000632; %波长(mm)
z0=1000;
L=N*pix; %CCD宽度(mm)
In(1:N,1:N)=f0(1:N,1:N);
%-----------------------------1-FFT重建开始
n=1:N;
x=-L/2+L/N*(n-1);
y=x;
[yy,xx] = meshgrid(y,x);
k=2*pi/h;
Fresnel=exp(i*k/2/z0*(xx.^2+yy.^2));
f2=In.*Fresnel;
Uf=fft2(f2,N,N);
Uf=fftshift(Uf);
L0=h*z0*N/L;
x=-L0/2+L0/N*(n-1);
y=x;
[yy,xx] = meshgrid(y,x);
phase=exp(i*k*z0)/(i*h*z0)*exp(i*k/2/z0*(xx.^2+yy.^2));
U0=Uf.*phase;%积分运算结果乘积分号前方相位因子
%-----------------------------1-FFT重建结束
If=U0*conj(U0);
Gmax=max(max(abs(U0)));
Gmin=min(min(abs(U0)));
figstr=strcat('重建物平面宽度=',num2str(L0),'mm');
figure(2),
imshow(abs(U0),[Gmin Gmax/1]),colormap(gray); xlabel(figstr);title('1-FFT物平面重建图像');
p=10;
while p
figure(3);
imshow(abs(U0),[Gmin Gmax/p]),colormap(gray);xlabel(figstr);title('1-FFT物平面重建图像');
p=input('Gmax/p,p=10?');
end;
Up=fft2(U0);
figure,imshow(abs(Up),[])
figure,plot(abs(Up(round(512)+1,:)));
hold on;
[1] 古德曼 ,傅里叶光学导论[M],北京:电子工业出版社,2006,第3版
[2]李俊昌,衍射计算及数字全息[M],北京:科学出版社,2014
[3]陈家壁,苏显渝.光学信息技术原理及应用[M].北京:高等教育出版社,2001