★ FEATURED ARTICLE 简介面向图像处理与计算成像方向的学习者这份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 31M 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*(x2*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((x35).^2 (y-25).^2 14^2) 0.85; % 圆 O((x-40).^2 (y20).^2 10^2) 0.6; % 小圆 O O .* exp(1i * 0.0015 * (x 2*y)); % 缓慢相位斜面 O_f F(O);这段代码的逻辑是先造振幅再乘相位最后做中心化傅里叶变换。F和IF必须成对使用F先ifftshift再fft2再fftshiftIF先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_x1) || any(idx_xN) || any(idx_y1) || any(idx_yN) 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要同奇偶。本文参数下N257、M63pad97正好居中。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(pad1:padM, pad1:padM) low0_f; O_est IF(O_f_est);这段代码先找中心 LED 索引led_x.^2 led_y.^2最小的那颗就是中心。low0是中心强度开方low0_f是它的频谱pad是本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
