光学成像中的MATLAB仿真：从衍射传播到PSF/MTF全解析

简介光学成像的MATLAB仿真代码包面向光学工程、计算成像方向的学生与科研人员用于理解光路建模、光线追迹、傅里叶变换、像差分析等核心概念涵盖从几何光学到傅里叶光学的典型实现。压缩包共5个文件整体仅133KB包含3个MATLAB脚本.m、1个说明文本.txt和1张细胞示例图片.jpg脚本可直接运行或在此基础上调整参数txt说明文档用于快速理解代码结构与运行流程。内容预览可看到细胞图像的成像仿真示例配合脚本便于对照学习直观呈现光线传播与成像结果也适合验证点扩散函数、调制传递函数等成像质量指标。目前已有24人学习适合需要快速上手光学成像仿真、验证理论算法的初学者可借助现成代码缩短搭建环境的时间并将注意力集中在成像原理与MATLAB实现技巧上。1. 光学成像的MATLAB仿真的真正起点是什么大多数第一次接触光学成像的MATLAB仿真的人会立即去找Image Processing Toolbox里的imfilter或fspecial把一张清晰图片和一个高斯核卷积然后当作“镜头退化”。这个流程能出图但对接下来的镜头设计、探测器采样、像差分析或复原算法验证几乎没有帮助因为你无法回答“为什么是这个模糊半径”“焦外为什么是这种形状”。光学成像的MATLAB仿真本质上是一个连续光场离散化和传播的问题用复振幅矩阵表示波前用衍射算子计算传播再用点扩散函数描述系统对点光源的响应。只要网格尺寸、波长、传播距离三者不满足采样条件FFT结果就会出现低频虚假结构甚至发散。这篇文章会按“衍射传播→PSF/MTF→整幅图像仿真”的顺序把每个环节的参数设定、代码和排错方法讲清楚。2. 建立标量衍射模型从菲涅尔传播到MATLAB实现2.1 为什么成像仿真默认采用标量衍射近似光学成像系统里的电磁场完整描述是矢量场但绝大多数镜头系统工作在近轴区域光束的偏振方向在传播中几乎不改变不同偏振分量之间的耦合可以忽略。此时把电场分量分别看作标量用复数振幅 (U(x,y)) 描述空间某点的场强和相位就能解释干涉、衍射、焦深和分辨率这些核心现象。这一近似在数值孔径NA小于0.7时足够准确如果仿真高NA显微物镜、超透镜或深亚波长光栅成像就需要转用RCWA或FDTD这类矢量方法。MATLAB里做标量衍射仿真第一步不是调函数而是确定光场的离散方式。一个二维光场 (U(x,y)) 被采样成 (N \times N) 复矩阵相邻采样点的实际间距是 (\Delta x)那么整个计算区域边长是 (N\Delta x)。成像仿真的很多错误都出在 (\Delta x) 与波长的量级关系上(\Delta x) 必须小于最小结构尺寸的一半同时还要保证传播过程中的相位变化不产生混叠。N 1024; % 网格点数用偶数是惯例 dx 10e-6; % 源平面采样间隔单位米 x (-N/2 : N/2-1) * dx; % 中心对称坐标 [X, Y] meshgrid(x); lambda 633e-9; % He-Ne激光波长 k 2*pi / lambda; % 波数逻辑说明这段代码建立了一个最基础的光场坐标系统之后所有的相位因子都基于X和Y计算。N取偶数是因为fftshift和ifftshift在偶尺寸下行为一致频率轴可以直接用-N/2 : N/2-1表示奇数会多出一个零频位置容易让人在坐标对齐上犯迷糊。dx的单位是米不要被图像处理的“像素”概念带偏后面计算目标平面坐标时dx直接进入衍射公式。参数说明内存占用是复矩阵每个元素16字节N2048时仅U就约64MB再加上X、Y两个坐标矩阵会到128MB以上。因此我一般建议先用1024试通流程确认坐标映射正确后再放大N否则调一次错一次时间都耗在等待上。2.2 菲涅尔单次FFT实现公式、代码与目标平面坐标标量衍射中最常用的算子是瑞利-索末菲衍射积分。在近轴条件下它可以简化为菲涅尔衍射积分。菲涅尔积分的离散实现有很多种最稳定也最容易写对的是“单次FFT形式”先给源平面乘一个二次相位因子再做一个FFT最后在输出平面乘另一个二次相位因子。输出平面坐标不是随便取的而是由空间频率和传播距离唯一决定(x_2 \lambda z f_x)。如果忽略这层坐标换算光斑位置和尺寸都会错得离谱。function [U2, dx2] fresnel_fft(U1, dx1, lambda, z) % 单次FFT菲涅尔衍射传播 % U1: 源平面复振幅矩阵 % dx1: 源平面采样间隔(m) % lambda: 波长(m), z: 传播距离(m) N size(U1, 1); L1 N * dx1; k 2*pi / lambda; x1 (-N/2 : N/2-1) * dx1; [X1, Y1] meshgrid(x1); % 源平面的二次相位因子 U1q U1 .* exp(1i * k / (2*z) * (X1.^2 Y1.^2)); % 中心化FFTifftshift把原点移到矩阵左上角fftshift再还原 Uf fftshift(fft2(ifftshift(U1q))); % 输出平面的空间频率坐标 fx (-N/2 : N/2-1) / L1; x2 lambda * z * fx; [X2, Y2] meshgrid(x2); % 输出平面的二次相位因子 U2 exp(1i*k*z) / (1i*lambda*z) ... .* exp(1i * k / (2*z) * (X2.^2 Y2.^2)) .* Uf; dx2 x2(2) - x2(1); % 输出平面采样间隔 end逻辑说明源平面二次相位因子 (e^{i\frac{k}{2z}(x_1^2y_1^2)}) 是菲涅尔积分的核心它表示球面波在到达观测点时的相位差。ifftshift和fftshift的配对一定不能乱写原始矩阵的中心在数组中央但FFT的零频在左上角所以要先ifftshift再fft2频谱算完后又用fftshift把零频移回中央这样后面的坐标轴中心对齐。输出平面的采样间隔 ( \Delta x_2 \lambda z / (N \Delta x_1)) 会自动变化z 越大(\Delta x_2) 越大也就是观察窗口尺寸 (N\Delta x_2) 越大。参数说明使用这个函数时dx1、lambda、z、N不是独立变量。单次FFT菲涅尔法的输出窗口尺寸和源窗口尺寸之间有一个采样约束。如果输出窗口小于你关心的物理区域说明要么减小dx1要么增大N要么改用角度谱法。当输出出现NaN或峰值出现在角落时先检查z是否太小再检查二次相位因子中是否有除以零。2.3 角度谱法近距离传播更稳健的替代方案当传播距离只有几毫米甚至几个波长时菲涅尔单次FFT的输出采样间隔会变得极小输出窗口缩小很多仿真结果往往是一团模糊。此时应该使用角谱法把源平面分解成平面波分量每个分量乘以一个相位延迟 (e^{i k_z z})再合成输出平面。角谱法的最大好处是输出平面的采样间隔和输入完全一样计算区域大小不变。function U2 angular_spectrum(U1, dx, lambda, z) % 角谱法传播保持采样间隔不变 % dx: 源平面和目标平面共同的采样间隔 N size(U1, 1); L N * dx; k 2*pi / lambda; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx); % 平面波传播因子的频域表达 H exp(1i * k * z .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 倏逝波在视野外直接置零避免sqrt负数造成NaN H((lambda*FX).^2 (lambda*FY).^2 1) 0; Uf fftshift(fft2(ifftshift(U1))); U2 ifftshift(ifft2(fftshift(Uf .* H))); end逻辑说明H就是频域传递函数每一项对应一个平面波方向。当空间频率平方和超过 (1/\lambda^2) 时对应的模式是倏逝波在传播距离远大于波长时衰减到可以忽略直接置零比保留NaN好得多。角谱法对距离没有下限限制但距离增大时高频振荡的传递函数会使数值误差累积所以它更适合近场和毫米级距离的传播。实际工程里我一般用下面的经验法则当 (z N\Delta x^2 / \lambda) 时用菲涅尔单次FFT当 (z 10 \cdot N\Delta x^2 / \lambda) 时用角度谱验证避免采样发散。两种算子得到的峰值位置和能量偏差在5%以内就说明网格参数没选错。传播算子输出采样间隔主要限制典型适用场景菲涅尔单次FFT(\lambda z / (N\Delta x_1))输出窗口随z变化远场、聚焦系统、夫琅禾费区角度谱法与输入相同 (\Delta x)大z下数值误差累积近场衍射、波导端面、短距离传播3. 仿真点扩散函数与MTF评价成像质量的核心3.1 从出瞳函数到PSF傅里叶光学的一行关系成像系统对点光源的响应是点扩散函数PSF它是成像质量的“指纹”。在非相干成像下强度PSF等于出瞳函数的傅里叶变换的模平方先做出瞳平面上的复振幅分布包含孔径形状和波前像差一次FFT得到焦平面电场振幅再取模平方就得到PSF。很多人用高斯函数拼一个PSF那只适合做图像处理演示不能用于设计镜头。出瞳函数 (P(x,y)) 通常包含两个部分振幅分布入瞳的形状圆形、环形、带遮挡和相位分布离焦、像散、球差等波前像差。MATLAB里把这两部分写进同一个复矩阵是PSF仿真的关键。function [psf, dx_psf] generate_psf(N, dx_pupil, lambda, focal, zernike) % 生成非相干成像系统的强度PSF % zernike: 结构体可包含 defocus, astig 等像差系数(以波长为单位) x (-N/2 : N/2-1) * dx_pupil; [X, Y] meshgrid(x); R_pupil dx_pupil * N/2; % 入瞳半径 rho sqrt(X.^2 Y.^2) / R_pupil; theta atan2(Y, X); pupil double(rho 1); % 圆形孔径 phase zeros(size(X)); if isfield(zernike, defocus) % Z4离焦: sqrt(3)*(2*rho^2-1) phase phase zernike.defocus * sqrt(3) * (2*rho.^2 - 1); end if isfield(zernike, astig) % Z5像散: sqrt(6)*rho^2*cos(2*theta) phase phase zernike.astig * sqrt(6) .* rho.^2 .* cos(2*theta); end phase 2*pi * phase .* pupil; % 单位转换为弧度 amp pupil .* exp(1i * phase); % 出瞳复振幅 E_focal fftshift(fft2(ifftshift(amp))); % 焦平面电场 dx_psf lambda * focal / (N * dx_pupil); % 焦平面采样间隔 psf abs(E_focal).^2; psf psf / sum(psf(:)); % 能量归一化 end逻辑说明像差系数以“波长”为单位比如zernike.defocus 0.25表示离焦波前差峰值为四分之一波长。相位项乘以2*pi是因为一个波长的光程差对应 (2\pi) 相位变化。fftshift(fft2(ifftshift(amp)))的移位方式和前面菲涅尔传播一致保证焦平面中心对应阵列中心。dx_psf的物理含义是焦平面上每个像素对应的真实尺寸它由波长和焦距决定(\lambda f / (N \Delta x_{pupil}))与探测器像素尺寸没有直接关系。调用示例zernike.defocus 0.25; % 1/4波长离焦 [psf, dx_psf] generate_psf(512, 4e-6, 633e-9, 5e-3, zernike); figure; imagesc(abs(psf(241:272, 241:272))); % 看中心区域 axis image; colormap hot;这个示例中入瞳半径是 (512/2 \times 4,\mu m \approx 1.024,mm)焦距5mm所以焦平面采样间隔大约0.6微米。如果你想和探测器像素尺寸3.45微米对应还需要对PSF做重采样这一步留给第4章。3.2 从PSF到MTFFFT的正确姿势与频率坐标换算MTF是光学系统对比度传递能力的频域描述它是PSF的傅里叶变换的模。很多人在这一步算出的MTF曲线中心不在1或者在原点处有相位跳变原因多半是忘了归一化或fftshift使用错误。正确做法是先对PSF做ifftshift再FFT把直流分量移到中心再取模。otf fft2(ifftshift(psf)); otf otf / sum(psf(:)); % 直流分量为1 mtf abs(otf); % MTF是模值 N size(psf, 1); fx (-N/2 : N/2-1) / (N * dx_psf); % 空间频率单位周期/米 figure; plot(fx(N/21:end), mtf(N/21, N/21:end)); % 沿x方向画MTF曲线 xlabel(空间频率 (cycles/mm)); ylabel(MTF); xlim([0, 1/(2*dx_psf)*1e-3]); % 奈奎斯特频率逻辑说明psf中心在矩阵中央ifftshift把中心搬到左上角FFT的结果零频也在左上角但显示时通常让零频在中央。用sum(psf(:))归一化是因为PSF的总能量已归一为1所以OTF在零频处自动为1如果PSF没有归一化这里必须补上。频率轴千万不能直接写(-N/2:N/2-1)/N那样得到的是归一化频率不是物理频率。要得到物理频率必须除以实际面阵尺寸 (N \Delta x_{psf})。3.3 不同像差对PSF和MTF的视觉影响只用一组像差系数看不到规律我把离焦、像散、球差三个典型结果放到同一张表里仿真时可以直接对照。像差类型Zernike项PSF特征MTF特征离焦(Z_4)中心能量分散出现圆环中频曲线提前下降整体截止频率降低像散(Z_5)水平和垂直方向不对称两个方向的MTF曲线分离容易出现“z”形交汇球差(Z_{11})中心亮核外带明显旁瓣中频区域出现非单调起伏不能简单用高斯拟合理解这个表的意义在于当你的光学成像的MATLAB仿真结果里PSF出现非高斯旁瓣时那不是噪声而是像差在起作用。比如模拟一个离焦0.4波长的镜头你会在PSF中心看到甜甜圈结构这是正常现象而不是程序写错。反过来如果图像退化仿真里用的PSF是完美高斯那说明你跳过了光学设计中最关键的一步。4. 仿真整个成像过程从场景到探测器图像4.1 PSF重采样让像素坐标和物理坐标对齐单个PSF算出来以后要和场景图像做卷积。但场景图像的单位是“像素”PSF的单位是“米”两者直接卷积没有物理意义。对于无穷远物体、有限口径透镜的成像系统物面一点在像面上的分布就是PSF且像面坐标已经由几何放大率决定。所以实际仿真时需要把PSF重采样到探测器像素尺寸上。dx_det 3.45e-6; % 探测器像元边长单位米 scale dx_psf / dx_det; % PSF采样间隔和像元尺寸的比值 % 重采样scale1表示PSF原始间隔大于像元向上插值scale1则向下抽稀 psf_r imresize(psf, scale, bilinear); psf_r psf_r / sum(psf_r(:)); % 重采样后重新归一化逻辑说明imresize缩放的是像素网格本质上相当于在物理坐标上做了线性重采样。这里用双线性足够因为PSF是平滑函数。如果dx_psf和dx_det相差超过5倍建议改用spline插值避免丢失PSF高频细节。重采样后必须再次归一化否则卷积后图像亮度会漂移。4.2 空间域卷积还是频域滤波按PSF尺寸决定重采样后的PSF如果小于64×64像素直接conv2是最快、最不容易出错的方式如果PSF接近甚至超过图像尺寸就应该用频域乘法。频域乘法的坑是循环卷积会造成边缘“卷绕”需要根据场景边缘是否重要来决定补零方式。function img_blur fft_convolve(img, psf) % 频域卷积等价于循环卷积使用时注意边界 M size(img, 1); N size(img, 2); H fft2(ifftshift(psf), M, N); % 将PSF填充到图像尺寸并FFT img_blur ifft2(fft2(img) .* H, symmetric); end逻辑说明fft2(ifftshift(psf), M, N)先把psf中心移到左上角然后自动补零或截断到M×N。symmetric告诉MATLAB输出应该接近实值把数值虚部直接清零避免imagesc显示时出现浮点尖刺。这个函数实现的是循环卷积图像左边缘会卷积到右边缘对于模拟光学成像我一般会先把场景边缘做一圈镜像扩展卷积完再裁剪得到无人工边缘的结果。空间域和频域的选择可以总结为下表方法适用PSF尺寸边界处理推荐场景conv2(im, psf, same)小于64×64零填充但边缘偏暗算法验证、小图仿真imfilter(im, psf, replicate)小于300×300边缘复制视觉自然自然图像退化仿真fft_convolve超过300×300必须手动补零或镜像扩展PSF很大、需要高速处理4.3 一个可运行的成像仿真主流程把上面的PSF生成、重采样、卷积、噪声、探测器采样放在一起就是一个比较完整的成像仿真链路。下面的脚本可以照抄。% 成像仿真主流程 scene im2double(imread(cameraman.tif)); scene scene / max(scene(:)); % 归一化到0-1 zernike.defocus 0.1; % 少量离焦 [psf, dx_psf] generate_psf(256, 8e-6, 550e-9, 5e-3, zernike); dx_det 3.45e-6; psf imresize(psf, dx_psf / dx_det, bilinear); psf psf / sum(psf(:)); % 卷积模拟光学模糊 img_blur conv2(scene, psf, same); % 光子噪声泊松分布scale越大噪声越低 scale 1e5; img_noise poissrnd(img_blur * scale) / scale; % 2x2像素合并模拟大像元 [M, N] size(img_noise); M M - mod(M, 2); N N - mod(N, 2); img_noise img_noise(1:M, 1:N); img_binned squeeze(mean(reshape(img_noise, [2, M/2, 2, N/2]), [1 3]));逻辑说明generate_psf生成的PSF默认在光学焦面上dx_det是探测器像元物理尺寸两者必须匹配。conv2(...,same)在PSF尺寸小时可以接受如果PSF的真实尺寸超过30像素建议换成4.2的频域实现。泊松噪声直接作用在光强上scale越大每个像素平均光子数越多噪声相对越小。像素合并模拟了物理探测器的分辨能力2×2合并后人眼看起来“更干净”但高频细节已经损失。参数说明这个脚本里的zernike.defocus 0.1对应的是一小段离焦量想要模拟严重失焦可以提高到0.5此时图像中心对比度会明显下降。受噪声影响后MTF的实际可分辨频率接近奈奎斯特频率的一半这属于正常现象。4.4 加快大场景仿真的两个方向GPU阵列和分块卷积4K、8K图像做光学成像仿真时conv2会非常慢。常见做法是先把PSF和场景都转成gpuArray用imfilter在GPU上跑如果GPU显存不够就使用blockproc分块。psf_g gpuArray(psf); scene_g gpuArray(scene); img_blur_g imfilter(scene_g, psf_g, replicate, same); img_blur gather(img_blur_g);逻辑说明imfilter在GPU上对浮点矩阵的速度通常比CPU快一个数量级但输入图像和PSF必须同时放进显存。4K灰度图大约32MB加上PSF和临时变量一般8GB显存够用。注意imfilter的replicate边界模式比零填充更接近真实传感器边界适合图像边缘对比度较高的场景。如果GPU内存紧张用blockproc按512×512分块block_size 512; border floor([size(psf,1) size(psf,2)]/2); img_blur blockproc(scene, [block_size block_size], ... (b) imfilter(b.data, psf, replicate), ... BorderSize, border, ... TrimBorder, true, ... PadPartialBlocks, true);逻辑说明BorderSize必须是PSF半宽以上否则块与块边界会出现一条明显的断层PadPartialBlocks用于处理图像宽高不能被块大小整除的情况。分块卷积的开销主要在边界重复计算块越大效率越高但单块太大时内存又吃紧一般512左右是平衡点。5. 用MTF/PSF验证仿真参数三个容易翻车的检查点5.1 能量守恒检查发散先从总功率开始查光学成像的MATLAB仿真出现发散时很多人直接怀疑算法但最常见原因是网格采样不够。最简单的检查是算传播前后的总能量能量变化超过5%说明频谱被截断或算子不匹配。E0 sum(abs(U1(:)).^2) * dx1^2; E1 sum(abs(U2(:)).^2) * dx2^2; fprintf(能量变化 %.2f%%\n, abs(E1-E0)/E0*100);能量变化大时逐级缩小dx并增大N直到变化稳定在1%以内。别指望能量完全不差FFT本身的离散误差会带来1%左右的波动。5.2 PSF重采样检查像素坐标和物理坐标是否对齐很多人在第4章把dx_psf和dx_det搞混导致图像模糊程度忽大忽小。验证方法是在重采样前后分别打印两个采样间隔确认重采样后的PSF半高宽在探测器上占了多少像素。fprintf(dx_psf%.3fum, dx_det%.3fum\n, dx_psf*1e6, dx_det*1e6);如果dx_psf远小于dx_det说明PSF在探测器上已经欠采样图像看起来像没退化如果远大于图像会过度模糊。正确的做法是让dx_psf和dx_det接近再通过imresize微调。5.3 相位混叠检查仿真发散的最常见来源二次相位因子在网格边缘变化很快如果相邻采样点之间的相位差超过(\pi)FFT会把高频折叠到低频形成错误的干涉条纹。这个检查应该在每次修改dx、z或lambda后执行。phase1 angle(U1q(1:end-1, :) .* conj(U1q(2:end, :))); if max(abs(phase1(:))) pi warning(相邻相位差超过pi请减小dx或增大N); end逻辑说明U1q是菲涅尔传播中源平面加过二次相位因子的复矩阵相邻两点相乘取相位差得到的就是该方向上的相位增量。如果任何一个方向超过(\pi)采样不满足奈奎斯特定理FFT会引入虚假低频成分。这个检查放在菲涅尔函数内部而不是放在主循环里能省下大半“仿真发散”的调试时间。本文还有配套的精品资源点击获取