☰
MATLAB傅里叶叠层成像仿真:低NA强度图重建高分辨率复振幅
2026/10/26 19:39:30 网站建设 项目流程

简介:面向图像处理与计算成像方向的学习者,这份MATLAB仿真材料围绕傅里叶叠层成像(FPM)的图像重建流程展开,适合用于理解多角度照明下的相位恢复、频域拼接与高分辨率重建等核心概念,也可作为课程设计或算法入门练习的参考。压缩包共1个文件,主体为.m脚本,大小约5KB,结构精简,便于直接阅读、运行和调试。已有2604人学习,说明其在FPM入门实践中具有一定参考价值。内容涉及模拟含振幅与相位信息的低分辨率图像、利用二维离散傅里叶变换获取部分傅里叶谱,并通过迭代优化更新高分辨率估计,配合图像显示观察中间结果与最终重建效果;读者可借此梳理FPM从测量建模、频域融合到收敛判断的实现思路,加深对傅里叶变换在成像系统中作用的理解。

1. 一组低 NA 强度图,为什么值得用 MATLAB 拆一遍

一组低 NA 物镜拍下来的 LED 阵列强度图,肉眼只能看到模糊的强度分布,真正的高频相位藏在不同照明角度对应的频域平移里。这套 MATLAB 傅里叶叠层图像重建仿真,先用虚拟高分辨率复振幅目标生成 15×15 的强度图栈,再用交替投影把低分辨率强度图拼回高分辨率复振幅,频谱平移、孔径低通、强度替换和迭代步长都能单独改。适合手里有算法验证任务、又不想从零搭光学平台的人,也适合课程设计里要把相位恢复过程跑通、还要能解释每个参数为什么这么设的场景。它不是只给一张漂亮结果图,而是把“强度图栈怎么来、频谱怎么平移、迭代为什么收敛”拆成可改参数的 MATLAB 脚本。

2. 傅里叶叠层的物理链路:LED 角度、频域平移与低 NA 采集

傅里叶叠层成像的核心链路并不复杂:LED 阵列中每颗 LED 以不同角度照明样品,样品出射场在物镜后焦面发生频谱平移,低 NA 物镜只让孔径内的低频分量通过,相机记录的是对应角度下的低分辨率强度图。换一颗 LED,就相当于把样品频谱平移到孔径位置再截取一次。不同角度下的强度图在频域里互相重叠,重叠区域提供冗余约束,交替投影利用这些约束把低 NA 孔径外的频谱补出来,最终得到高分辨率复振幅。MATLAB 仿真的价值在于每个物理量都能落到像素:LED 坐标、照明角度、频谱平移量、孔径半径、强度图尺寸、重叠率,全都可以单独改。

2.1 为什么用频谱平移,而不是空域插值超分

空域插值超分只利用单张图的像素相关性,缺的是高频相位信息。傅里叶叠层不靠插值猜细节,它把不同照明角度下的低分辨率强度图看成同一个高分辨率频谱在不同平移位置的低通采样。每张强度图对应一个频域子区域,子区域之间有重叠,重叠率越高,频谱拼接的约束越强。常见做法是让相邻 LED 的频谱平移量小于孔径直径的一半,重叠率保持在 60% 到 80% 之间。重叠太少,频谱拼图会出现缝隙;重叠太多,采集帧数变多,迭代时间上升,但对噪声更稳。MATLAB 仿真里可以先把重叠率算出来,再决定 LED 间距和阵列高度。

2.2 参数落到像素:shift 与 aperture 的公式

在仿真里,LED 角度用tanθ_x = x_led / h近似,小角度下sinθ_x ≈ x_led / h。空间频率平移量是Δf_x = sinθ_x / λ,离散频域采样间隔是df_x = 1/(N dx_obj),所以频谱平移像素数为s_x = x_led N dx_obj / (h λ)。孔径半径对应物镜截止频率,r_ap = NA N dx_obj / λ。这两个公式决定了重建画布、强度图尺寸和最大 LED 角度。下面这段 MATLAB 脚本把参数表直接跑出来,顺便检查边缘 LED 会不会超出画布。

% 傅里叶叠层仿真参数:先算频谱平移和孔径半径 lambda = 0.532e-6; % 波长,m NA = 0.13; % 低 NA 物镜数值孔径 N = 257; % 重建画布像素数,取奇数便于 fftshift dx_obj = 0.5e-6; % 物面采样间隔,m led_pitch = 4e-3; % LED 间距,m led_height= 80e-3; % LED 平面到样品距离,m led_num = 15; % 单边 LED 数,15x15 r_ap = round((NA/lambda) * N * dx_obj); % 孔径半径,频域像素 M = 2*r_ap + 1; % 单张低分辨率强度图尺寸 [led_x, led_y] = meshgrid(((1:led_num)-ceil(led_num/2))*led_pitch); led_x = led_x(:); led_y = led_y(:); shift_x = round((led_x / led_height) * N * dx_obj / lambda); shift_y = round((led_y / led_height) * N * dx_obj / lambda); max_shift = max(abs([shift_x; shift_y])); fprintf('r_ap=%d, M=%d, max_shift=%d, 画布半宽=%d\n', ... r_ap, M, max_shift, floor(N/2)); if max_shift + r_ap > floor(N/2) warning('最大平移加孔径超过画布半宽,边缘 LED 会被截断'); end

这段代码的逻辑很直接:r_ap是低 NA 孔径在频域里的半径,M是相机看到的低分辨率图尺寸,shift_x和shift_y是每颗 LED 对应的频域平移像素。led_x和led_y用meshgrid生成中心对称的 15×15 坐标,中心 LED 坐标为 0。最后的边界检查很重要,max_shift + r_ap超过floor(N/2)时,边缘 LED 的频谱区域会被画布截断,重建结果会在高频出现暗角。参数改动时,led_height调大,角度变小,shift变小;NA调大,r_ap变大,低分辨率图尺寸变大;N调大,频域采样变密,shift和r_ap同步变大。常见做法是先固定N和dx_obj,再根据物镜 NA 算r_ap,最后用led_pitch和led_height调重叠率。

参数数值作用改动后的连锁反应
lambda532 nm照明波长改小后频谱平移和孔径半径都变大
NA0.13物镜数值孔径调大后r_ap变大,低分辨率图变大
N257重建画布像素数调大后频域采样变密,计算量上升
dx_obj0.5 μm物面采样间隔调小后视场变小,shift和r_ap变小
led_pitch4 mmLED 间距调大后相邻频谱平移变大,重叠率下降
led_height80 mmLED 到样品距离调大后角度变小,shift变小
led_num15单边 LED 数调大后帧数增加,最大角度变大

按上表算出来,r_ap = 31,M = 63,相邻 LED 频谱平移约 12 像素,孔径直径 62 像素,重叠率约 80%,最大平移约 84 像素,边缘余量够用。这个参数组合适合先跑通全流程:15×15 共 225 张 63×63 强度图,双精度存储约 7 MB,普通笔记本就能跑。调试阶段可以先把led_num改成 5 或 7,用中心区域 LED 跑小阵列,确认频谱平移方向和重建逻辑没问题,再上全阵列。这样做的原因是全阵列迭代一次要遍历 225 张图,参数错了很难定位是方向问题还是归一化问题。

2.3 低 NA 采集在仿真里怎么表达

低 NA 采集不是简单把高分辨率图做高斯模糊。正确的表达是:高分辨率物体频谱平移到对应 LED 角度后,只保留孔径内的频谱,再逆变换得到低分辨率复振幅,相机取强度。也就是说,每张强度图都是高分辨率频谱的一个带限子区域。仿真里不需要真的模拟透镜传播,只要在频域做平移和截取,就能得到和低 NA 物镜一致的带限观测。这样做的好处是可控:孔径函数、平移量、噪声、归一化都能单独改。常见做法是把高分辨率复振幅目标放在N×N画布中心,频谱用fftshift(fft2(ifftshift(O)))中心化,然后对每颗 LED 取M×M的频谱子块。子块的中心偏移就是shift_x和shift_y,子块大小由r_ap决定。这个表达方式也让重建端和仿真端共用同一套索引,减少坐标不一致导致的翻车。

3. MATLAB 仿真骨架:从目标复振幅到 15×15 强度图栈

仿真骨架分三步:生成带相位的高分辨率复振幅目标,按 LED 角度生成低分辨率强度图栈,保存真值和参数。目标不能只用二值图像,最好带缓慢相位斜面和不同反射率区域,这样重建出来的相位才有东西可看。强度图栈要统一归一化,不能每张图单独归一化到峰值 1,否则不同 LED 之间的相对亮度丢失,重建时幅度约束会变弱。保存时把I_stack、O、led_x、led_y、shift_x、shift_y、r_ap、M、N、dx_obj一起存进 MAT 文件,后面调参和复现都靠这份文件。

3.1 生成带相位的高分辨率目标

目标复振幅用几何图形拼出来,横条、圆、小圆分别给不同振幅,再乘一个缓慢相位斜面。相位斜面不要太大,控制在0.0015*(x+2*y)这种量级,避免相位包裹太快导致逆变换出现剧烈条纹。目标放在画布中心,坐标用(0:N-1)-floor(N/2),这样中心像素对应零坐标。频谱用F和IF两个匿名函数封装,所有正变换和逆变换都走这两个入口,后面排查fftshift与ifftshift混用问题时省事。

F = @(x) fftshift(fft2(ifftshift(x))); IF = @(X) fftshift(ifft2(ifftshift(X))); N = 257; dx_obj = 0.5e-6; [x, y] = meshgrid((0:N-1)-floor(N/2)); O = zeros(N); O(abs(x)<18 & abs(y)<6) = 1.0; % 横条 O((x+35).^2 + (y-25).^2 < 14^2) = 0.85; % 圆 O((x-40).^2 + (y+20).^2 < 10^2) = 0.6; % 小圆 O = O .* exp(1i * 0.0015 * (x + 2*y)); % 缓慢相位斜面 O_f = F(O);

这段代码的逻辑是先造振幅,再乘相位,最后做中心化傅里叶变换。F和IF必须成对使用:F先ifftshift再fft2再fftshift,IF先ifftshift再ifft2再fftshift,这样频谱零频始终在中心。参数N取奇数 257,是为了让中心像素唯一,fftshift和ifftshift在奇数尺寸下不会出现半像素错位。dx_obj只参与频谱平移和孔径半径计算,不直接进入fft2,但后面算shift_x和r_ap时要用。目标振幅范围控制在 0 到 1,相位范围控制在几弧度以内,重建结果不会因为相位跳变太剧烈而出现条纹伪影。

3.2 生成每颗 LED 的低分辨率强度图

生成强度图时,先按第 2 章算好的shift_x、shift_y和r_ap,从O_f里取M×M子块。子块中心是画布中心加平移量,索引范围是cx + shift_x(j) + (-r_ap:r_ap)。取完子块后做逆变换得到低分辨率复振幅,取abs().^2得到强度。每张强度图先存原始值,循环结束后统一除以全局最大值,保证相对亮度一致。如果仿真里想模拟不同 LED 亮度差异,可以在循环里乘一个led_gain(j),但重建端要知道这个增益,否则幅度约束会偏。

r_ap = 31; M = 2*r_ap + 1; I_stack = zeros(M, M, numel(led_x)); for j = 1:numel(led_x) cx = floor(N/2) + 1; cy = cx; idx_x = cx + shift_x(j) + (-r_ap:r_ap); idx_y = cy + shift_y(j) + (-r_ap:r_ap); if any(idx_x<1) || any(idx_x>N) || any(idx_y<1) || any(idx_y>N) error('第 %d 个 LED 的频谱区域越界', j); end patch = O_f(idx_y, idx_x); low = IF(patch); I_stack(:,:,j) = abs(low).^2; end I_stack = I_stack / max(I_stack(:)); save('fpm_stack.mat','I_stack','O','led_x','led_y', ... 'shift_x','shift_y','r_ap','M','N','dx_obj');

代码里的cx和cy是画布中心索引,idx_x和idx_y是当前 LED 的频谱子块范围。patch是M×M的频域子块,IF(patch)得到M×M的低分辨率复振幅,abs(low).^2是相机强度。边界检查放在取子块之前,越界直接报错,不要用circshift自动环绕,环绕会把另一侧的高频混进来,重建结果出现周期性拼接缝。save把真值O和所有参数一起保存,后面算相位误差和做消融都直接用这份文件。参数说明:r_ap决定低分辨率图尺寸,M必须是奇数,I_stack第三维是 LED 数量,led_x和led_y的顺序要和shift_x、shift_y严格对应。

3.3 强度图栈的检查清单

生成完不要急着跑重建,先做三个检查。第一,看中心 LED 强度图是不是目标低频模糊像,如果中心图都是黑的,说明O_f子块索引偏了。第二,看边缘 LED 强度图有没有明显变暗或出现硬截断,如果边缘图一侧发黑,说明max_shift + r_ap接近画布半宽,需要增大N或减小led_height。第三,看I_stack的全局最大值和最小值,最小值接近 0 是正常的,但如果大量像素为 0,说明归一化或目标振幅有问题。常见做法是抽三张图显示:中心 LED、最左 LED、最右 LED,确认强度分布随角度平移连续变化。这样做的原因是重建端靠强度图栈提供幅度约束,强度图栈本身有截断或归一化错误,后面调迭代步长救不回来。

4. 交替投影重建:频谱约束、强度替换与步长参数

重建端的输入只有强度图栈和 LED 平移量,目标复振幅O只用来做验证,不能进迭代。初始化用中心 LED 的低分辨率强度补零到N×N频域,得到初始高分辨率复振幅。每轮迭代随机打乱 LED 顺序,对每颗 LED 做四件事:从当前频谱取对应子块,逆变换得到低分辨率复振幅,用测量强度替换振幅并保留相位,正变换后把差值按松弛因子更新回频谱子块。这个过程就是交替投影:空域幅度约束和频域孔径约束来回投影,重叠区域的冗余信息把孔径外频谱逐步补出来。

4.1 初始化:为什么用中心 LED 补零

中心 LED 照明角度接近零,频谱平移量最小,它的低分辨率强度图对应样品频谱的低频部分。把这张图的振幅开方后做傅里叶变换,得到M×M的低频频谱,再补零到N×N画布中心,就得到初始高分辨率频谱。补零不是简单插值,它把低分辨率频谱放到高频画布的中心,高频区域初始为零,后面靠其他 LED 的平移子块逐步填充。初始相位可以设为零,也可以加随机相位,但随机相位会让前几轮残差震荡,调试阶段建议零相位起步。补零位置必须居中,pad = (N-M)/2必须是整数,所以N和M要同奇偶。本文参数下N=257、M=63,pad=97,正好居中。

load('fpm_stack.mat'); F = @(x) fftshift(fft2(ifftshift(x))); IF = @(X) fftshift(ifft2(ifftshift(X))); [~, cidx] = min(led_x.^2 + led_y.^2); % 找中心 LED low0 = sqrt(I_stack(:,:,cidx)); low0_f = F(low0); pad = (N - M)/2; O_f_est = zeros(N); O_f_est(pad+1:pad+M, pad+1:pad+M) = low0_f; O_est = IF(O_f_est);

这段代码先找中心 LED 索引,led_x.^2 + led_y.^2最小的那颗就是中心。low0是中心强度开方,low0_f是它的频谱,pad是

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询