Zernike相称显微镜的原理与matlab仿真

fightandstrive

已于 2024-01-23 11:54:14 修改

阅读量2.1k

点赞数 8

分类专栏：傅里叶光学&MATLAB编程文章标签：光学成像傅立叶分析 matlab 信号处理

于 2023-02-23 22:16:38 首次发布

本文链接：https://blog.csdn.net/qq_44791407/article/details/129188070

版权

傅里叶光学&MATLAB编程专栏收录该内容

6 篇文章 78 订阅

订阅专栏

1. Zernike相称显微镜的由来

1.1 对透明样本成像困难

人类生活在自然界，对世界充满好奇心，不断的探索使得人类能够逐渐掌握自然界的规律，并认识到生命的美丽。得以于16世纪显微镜的发明，人类也认识到生命起源于微小的生物大分子和细胞。然而，想要欣赏到微观世界的美丽并不容易，大部分的细胞和组织都是透明的，对光的吸收很弱，成像的对比度（最大灰度与最小灰度值之差）很低，这一难题困扰了生物学家和材料学家诸多年。虽然依靠化学试剂染色和荧光蛋白标记可以提高成像对比度，并且能够观察到彩色效果，但是，细胞也会面临被杀死的风险。

1.2 相位成像

人眼或者说相机能够获取包含样本信息的图像，是因为样本在被光照射之后，对光有所吸收且吸收有差异，使得到达人眼视网膜的光振幅不同。如果吸收没有差异，那么人眼看到的要么是一片白，要么一片灰，要么一片黑。

当光照射透明样本时，振幅的变化却相当小，人眼很难观察到样本信息。这时，能够利用的只有光波的相位，但是，无论是人眼还是现有的光电传感器，都很难直接观察到光波几百Thz频率的相位信息。这时，通过将光波的相位信息变为振幅信息的相位成像方法成为了救世之主。而Zernike相称就是一种非常优秀的相位成像方法，1953年的诺贝尔奖足以说明它的影响力。

2. Zernike相称显微镜的原理

2.1 对透明样本成像困难的物理解释

假定利用显微镜观察一个透明样本，照明光波为平面光波，样本的复振幅透过率为

$g(x,y)=exp[i\varphi (x,y)]$ (1)

其中， $\varphi (x,y)$ 是样本的相位， $(x,y)$ 是样本平面的二维空间坐标， $i$ 是单位虚数。因为透明样本吸收几乎没有，故振幅透过率为1（即e指数前系数为1），与此同时，因为透明样本在显微镜下观察时，通常被放置在与其折射率相匹配的培养液或者做成微米级别的切片，其相位一般很小，故可以对 $g(x,y)$ 作泰勒展开，得到

$g(x,y)=1+i\varphi (x,y)-\frac{1}{2}\varphi (x,y)^{2}-\frac{1}{6}i\varphi (x,y)^{3}+...\newline =1+i\varphi (x,y)$ (2)

使用普通显微镜直接观察时，为简单起见，假定方法倍数为1，不存在衍射和像差，那么像的强度分布为

$I(x',y')=\left | 1+i\varphi (x',y')\right |^{2}=1+\varphi (x',y')^{2}\approx 1$ (3)

其中， $(x',y')$ 时像平面的二维空间坐标。由此可见，拍摄到的图像对比度为0，要么是一片白，要么一片灰。

2.2 Zernike相称的物理解释

1935年，Zernike提出了观察位相物体的相衬法。Zernike认为，之所以观察不到相位引起的光强变化，是因为直接透过样本的强背景光（在公式(2)中实数1）和相位的弱衍射光 $\varphi (x,y)$ 之间存在 $\pi/2$ 的相位差，所以没办法干涉成像，那么如果他们能够同相位，就能观察到相位引起的光强变化。

但是，因为透过样本的背景光和衍射光几乎是耦合在一起，必须先想办法把它们分开。办法也很简单，加透镜即可，如图1所示，光源S发出光经过透镜L1准直之后照射到物平面P1上，直接透过的光会经过L2汇聚到焦点附近（零频背景光），而衍射光在焦平面的分布分散（低中高频衍射光），这时就非常适合在焦平面P2放置一个相位板，以改变零频和其他频率成分之间的相位关系。相位板的分布函数为

$H(f_{x},f_{y})=\left\{\begin{matrix} \pm j, f_{x}=f_{y} =0 \\ 1, other & \end{matrix}\right.$ (4)

2.3 正相称

考虑焦平面，光场分布为投射光波的频谱，即对式(2)作傅里叶变换，得到

$G(f_{x},f_{y})=\delta (f_{x},f_{y})+i\Phi (f_{x},f_{y})$ (5)

乘上相位板的分布函数，即对第一项附加 $\pi /2$ 的相位延迟，得到

$G(f_{x},f_{y})=i[\delta (f_{x},f_{y})+\Phi (f_{x},f_{y})]$ （6）

再作一次傅里叶变换，就能够得到像面上的光强分布为

$I(x',y')=\left | 1+\varphi (x',y')\right |^{2}=\\ 1+2\varphi (x',y')+\varphi (x',y')^{2}\approx 1+2\varphi (x',y')$ （7）

可以看到，此时像的强度与物体的相位成正比，称为正相称。

2.4 负相称

同样的，如果附加的相位是 $-\pi /2$ ，那么滤波后的频谱就变成了

$G(f_{x},f_{y})=i[-\delta (f_{x},f_{y})+\Phi (f_{x},f_{y})]$ （8）

再作一次傅里叶变换，就能够得到像面上的光强分布为

$I(x',y')=\left | -1+\varphi (x',y')\right |^{2}=\\ 1-2\varphi (x',y')+\varphi (x',y')^{2}\approx 1-2\varphi (x',y')$ (9)

此时，此时像的强度与物体的相位成反比，称为负相称。

2.5 中心暗场法

正相称和负相称都使得相位信息可视化，但是，强度变化 $2\varphi (x',y')$ 都是叠加在一个均匀背景‘1’上，此时的对比度为

$C=\frac{I_{max}(x',y')-I_{min}(x',y')}{I_{max}(x',y')+I_{min}(x',y')}\approx 2\varphi (x',y')$ (10)

可以发现，对比度仍然不高，想要提高对比度，就必须压缩零频。为解决这个问题，相位板除了对零频添加相位延迟之外，还添加了吸收项 $\alpha$ ，像的频谱就变为

$G(f_{x},f_{y})=\pm i[\alpha \delta (f_{x},f_{y})\pm \Phi (f_{x},f_{y})]$ （11）

再作一次傅里叶变换，就能够得到像面上的光强分布为

$I(x',y')=\left | \alpha \pm \varphi (x',y')\right |^{2}=\\ 1\pm 2\alpha \varphi (x',y')+\varphi (x',y')^{2}\approx \alpha ^{2} \pm 2\alpha\varphi (x',y')$ （12）

此时的对比度为 $2 \varphi (x',y')/ \alpha$ ，所以，想要提高对比度，只需要减小 $\alpha$ ，当 $\alpha$ 接近于0时，像的强度为

$I(x',y')= \varphi^{^{2}} (x',y')$ （13）

此时，称为中心暗场法，这也是Zernike相称显微镜所采用的方法。

2.6 思考

如果添加的相位是从0到 ${\color{Red} \pi}$ 逐渐变化，对比度会如何变化？不妨推导和仿真一下，很有趣的。

3. Zernike相称显微镜的仿真

3.1 代码

clc;
clear all;
close all;
%% (a) 初始化物体
jpg_iamge = imread('Rem.jpeg'); %读取相位图像，可以随便换，雷姆，狂三啥的
gray_iamge = double(imresize(im2gray(jpg_iamge),[250,250]));
gray_iamge = flipud(gray_iamge); %imagesc()函数的通病，需要反转一下，用imshow()就没事
gray_iamge = gray_iamge./max(max(gray_iamge));
[M,N]=size(gray_iamge); %get image sample size
phase=sqrt(gray_iamge)*pi/100;
ug=exp(1i*phase);
L=0.3e-3; %image plane side length (m)
dx=L/M; %sample interval (m)
x=-L/2:dx:L/2-dx; 
y=x;
%% (b) 普通成像
figure(1) %check ideal image
imagesc(x,y,abs(ug).^2);
colormap('gray'); colorbar;
xlabel('x (m)'); ylabel('y (m)');title('普通成像');
axis square
axis xy
figure(2) %horizontal image slice
vvalue=0; %select row (y value)
vindex=round(vvalue/dx+(M/2+1)); %convert row index
Ig2 = abs(ug).^2;
plot(x,Ig2(vindex,:)','r:');
xlabel('x (m)'); ylabel('Irradiance');
title('普通成像;(y=0m)');

%% (c) 正相称
Ug = fftshift(fft2(fftshift(ug)));
beta = pi/2;
Phase_modulation = exp(1i*beta);
Ug2 = Ug;
Ug2(M/2+1,M/2+1) = Ug2(M/2+1,M/2+1)*Phase_modulation;
uo = fftshift(ifft2(fftshift(Ug2)));
Io = abs(uo).^2;
figure(3)
imagesc(x,y,Io);
colormap('gray'); colorbar;
xlabel('x (m)'); ylabel('y (m)');title('正相称;\beta = \pi/2');
axis square
axis xy
vvalue=0; %select row (y value)
vindex=round(vvalue/dx+(M/2+1)); %convert row index
figure;
plot(x,Io(vindex,:)','r:');
xlabel('x (m)'); ylabel('Irradiance');
title('正相称剖线;\beta = \pi/2;y=0m');
%% (d) 负相称
beta = 3*pi/2;
Phase_modulation = exp(1i*beta);
Ug2 = Ug;
Ug2(M/2+1,M/2+1) = Ug2(M/2+1,M/2+1)*Phase_modulation;
uo = ifftshift(ifft2(Ug2));
Io = abs(uo).^2;
figure;
imagesc(x,y,Io);
colormap('gray'); colorbar;
xlabel('x (m)'); ylabel('y (m)');title('负相称;\beta = 3\pi/2');
axis square
axis xy
figure;
plot(x,Io(vindex,:)','r:');
xlabel('x (m)'); ylabel('Irradiance');
title('负相称剖线;\beta = 3\pi/2;y=0m');
%% (e) 中心暗场法
beta = pi/2;
Phase_modulation = 0.01*exp(1i*beta);
% Phase_modulation = 0*exp(1i*beta);
Ug2 = Ug;
Ug2(M/2+1,M/2+1) = Ug2(M/2+1,M/2+1)*Phase_modulation;
uo = ifftshift(ifft2(Ug2));
Io = abs(uo).^2;
figure
imagesc(x,y,Io);
colormap('gray'); colorbar;
xlabel('x (m)'); ylabel('y (m)');title('中心暗场法;\beta = 0.01exp(\pi/2)');
axis square
axis xy
figure;
plot(x,Io(vindex,:)','r:');
xlabel('x (m)'); ylabel('Irradiance');
title('中心暗场法;\beta = \beta = 0.01exp(\pi/2);y=0m');