简介基于Matlab的光场相机重聚焦算法演示方案适合正在学习计算摄影、光场成像或数字重聚焦的本科生、研究生及图像处理爱好者。方案围绕光场记录空间中光线方向与辐射强度这一核心思想演示了如何在完成拍摄后通过计算光学传输模型自由调整焦点从而突破传统光学系统单一焦深的局限。资源包共18个文件压缩后约13.59MB含3个.m源码文件、9张PNG结果对比图、md/txt说明文档、zbak备份文件及一个备用zip压缩包目录层次清晰。已有53人学习浏览。借助这套内容读者可以直接运行主程序和重聚焦函数观察不同焦点位置下的成像变化并对照结果图理解算法参数与效果之间的关系说明文档和备份文件也降低了复现门槛便于后续扩展成自己的光场成像实验项目。1. 光场相机拍下的不是一张照片而是一叠可以事后重聚焦的光场“拍的时候不对焦出片的时候再对焦”是光场相机最反直觉的地方。传统相机在按下快门前必须决定焦平面拍虚了只能靠后期去模糊硬拉。光场相机在主镜头和传感器之间塞了一层微透镜阵列入射光线的方向和位置被同时记录下来单次曝光得到的 RAW 图里实际上混叠着几十个视角的子图像。所谓重聚焦算法就是把这组子图像按不同深度重新投影并叠加让焦点在电脑上“滑”到任意深度平面。本文要做的就是把这套东西在 Matlab 里完整跑通光场数据怎么组织、两种主流重聚焦算法怎么实现、参数怎么定、演示界面怎么搭。适合做光场成像、三维重建方向的研究生也适合把图像处理课程设计做出区分度的本科生。2. 光场数据在 Matlab 里怎么摆微透镜图像、子孔径图像与光场重构2.1 从传感器 RAW 到 4D 光场微透镜阵列在做一次空间-角度复用光场用 4D 函数 L(u, v, s, t) 描述其中 (s, t) 是光线到达主镜头平面的位置(u, v) 是光线方向。微透镜阵列的作用是把这两个维度折叠到同一块二维传感器上每个微透镜覆盖一小块像素像素在微透镜内的坐标对应方向 (u, v)微透镜在阵列中的位置对应空间坐标 (s, t)。因此传感器上的一张图是 4D 光场的混叠投影不能直接当普通照片用。要让重聚焦算法能处理第一步是把微透镜图像“拆”成子孔径图像堆。子孔径图像相当于用一个很小的光圈从不同角度看同一个场景视角之间只存在视差这正是重聚焦算法需要的输入。拆的时候有两个关键尺寸微透镜间距每个微透镜覆盖几个像素常见 10~16 像素和微透镜阵列相对传感器的旋转角装配误差通常小于 1 度但必须补偿。数据处理方向上光场数据有三种常见的组织方式Matlab 里对应三种 array shape数据组织方式物理含义Matlab 数据形状典型用途微透镜图像堆传感器原始图按微透镜网格分块[pitch, pitch, M, N]白图像校正、RAW 解码子孔径图像堆固定 (u,v)按微透镜坐标取像素[nT, nS, nV, nU]重聚焦、视差估计EPI 极平面图固定 v 和 s切出 (u, t) 平面[nV, nT]深度估计、遮挡分析重聚焦算法全部建立在子孔径图像堆上所以接下来的重构代码以它为目标输出。2.2 在 Matlab 里把微透镜图像重构为子孔径图像堆真实相机的 RAW 解码涉及去马赛克、白图像校正、微透镜中心拟合环节很长。我这里给出工程上最常用的重构骨架拿到白平衡后的微透镜图 mli已知微透镜阵列网格数为 M×N、微透镜间距 pitch、旋转角 rotate按几何投影把每个微透镜下的像素块重排成子孔径图像。function lf MLICoord2Subview(mli, M, N, pitch, rotate) % MLICoord2Subview: 微透镜图像 - 子孔径图像堆 % mli: [h, w, 3] 白平衡后的微透镜图像 % M, N: 微透镜阵列的行、列数 % pitch: 每个微透镜覆盖的像素数直径 % rotate: 微透镜阵列相对传感器的旋转角单位度 % 输出 lf: [pitch, pitch, M, N, 3]单精度 [h, w, ~] size(mli); lf zeros(pitch, pitch, M, N, 3, single); cx (w 1) / 2; cy (h 1) / 2; for i 1:M for j 1:N % 1. 理想网格坐标 x0 (j - (N1)/2) * pitch cx; y0 (i - (M1)/2) * pitch cy; % 2. 补偿旋转绕传感器中心旋转 xr x0 * cosd(rotate) - y0 * sind(rotate); yr x0 * sind(rotate) y0 * cosd(rotate); x0 xr cx; y0 yr cy; % 3. 截取该微透镜覆盖的像素块 xs round(x0 - pitch/2 1 : x0 pitch/2); ys round(y0 - pitch/2 1 : y0 pitch/2); if any(xs 1) || any(ys 1) || any(xs w) || any(ys h) continue; % 越界透镜直接跳过真实数据边缘常见 end lf(:, :, i, j, :) mli(ys, xs, :); end end end这段代码的核心是坐标映射先算每个微透镜中心在理想网格上的位置再做一次绕中心的旋转补偿最后用 round 取整截块。rotate 参数在真实设备上通常由白图像white image的质心拟合得到如果没有标定数据可以先设 0 跑通流程再微调。输出 lf 的五个维度分别是微透镜内 y、微透镜内 x、微透镜行、微透镜列、RGB 通道。注意这里没有做子像素插值。pitch 是整数时误差可控如果 pitch 带小数比如 11.4应该用 interp2 做双线性采样否则重聚焦结果会出现固定条纹噪声。2.3 没有相机数据时用一张深度图生成仿真光场真实光场数据文件如 Lytro 的 .lfr需要专用解码库调试阶段不方便。我一般用一张纹理图和一张视差图先仿真出子孔径图像堆这样后面调重聚焦算法时每个参数的效果都是已知的。function lf SimulatedLF(textureMap, disparityMap, nU, nV) % SimulatedLF: 纹理图 视差图 - 子孔径图像堆 % textureMap: [h, w, 3] 纹理可来自imread % disparityMap: [h, w] 视差前景大、背景小范围建议 [0, 2] % nU, nV: 角度分辨率常用 5x5 或 7x7 % 输出 lf: [h, w, nV, nU, 3] [h, w, c] size(textureMap); lf zeros(h, w, nV, nU, c, single); [XX, YY] meshgrid(1:w, 1:h); for v 1:nV for u 1:nU % 归一化方向偏移范围约 [-1, 1] du (u - (nU1)/2) / nU; dv (v - (nV1)/2) / nV; % 视差越大偏移越大模拟不同深度在视角间的位移 Xs XX - du .* disparityMap .* 20; Ys YY - dv .* disparityMap .* 20; for ch 1:c lf(:, :, v, u, ch) interp2(XX, YY, ... double(textureMap(:, :, ch)), Xs, Ys, cubic, 0); end end end end生成原理是反向映射对每个视角 (u,v)把纹理图按视差图缩放后的坐标重采样。disparityMap 乘以 20 是把视差映射到像素位移的经验值实际场景里这个系数由标定基线决定。有了这个函数就可以用任意一张彩色图和它的语义分割深度图生成测试光场后续所有算法调试都不再依赖外部数据。3. 重聚焦算法为什么成立重投影模型与傅里叶切片定理3.1 数字重聚焦的几何直觉剪切、平移、再叠加光场重聚焦最常见的理论基础是“剪切-叠加”。假设场景中某个点深度为 d它在不同子孔径图里的像素位置差视差与子孔径坐标 (u,v) 成正比。要把焦点从当前传感器平面移动到深度为 α 的平面等价于把每个子孔径图按视差方向平移再累加求平均。用公式表达设参考焦点深度为 d_ref目标深度为 d重聚焦系数b 1 - d_ref / d那么第 (u,v) 个子孔径图的平移量是 b*(u-u_c, v-v_c)其中 (u_c,v_c) 是中心视角坐标。所有视角图平移后叠加得到的就是聚焦在深度 d 的重聚焦图像。b 为正表示焦点后移为负表示焦点前移。这个模型简单、直观也是实现在时域最直接的方式。时域平移加和的优点是物理意义清楚对遮挡和边缘的处理比较自然缺点是每个视角图都要做一次插值平移计算量与视角数成正比。对 7×7 视角、单张 512×512 的彩色图Matlab 循环实现大约要几百毫秒实时拖动会卡顿。3.2 傅里叶切片定理把重聚焦变成一次 2D 切片Ren Ng 对数字重聚焦的经典贡献是发现它和 X 射线断层扫描里的傅里叶切片定理同源一个 4D 光场的 2D 傅里叶切片对应某个深度的重聚焦图像。换句话说对光场 L(u,v,s,t) 做 4D 傅里叶变换然后在频域取一个经过原点的 2D 平面再做 2D 逆变换就得到重聚焦结果。频域切片的优势是复杂度和空间分辨率无关只取决于 FFT 和一次切片插值。对需要大量渲染焦栈几十个深度平面的应用频域法比时域平移加和快一个数量级。代价是频域处理对采样率敏感角度分辨率不足时切片会产生混叠表现为重聚焦图像上的周期性波纹。实际工程中也可以用“频率相位斜坡叠加”来等效实现避免 4D 数组切片的内存开销。3.3 两种路线的选型时域精度、频域速度、混合方案算法路线核心操作速度内存占用适用场景时域平移加和每视角 imtranslate 累加慢O(nU*nV)低逐视角处理少数量聚焦、验证算法、精度敏感傅里叶切片4D FFT 2D 切片 IFFT快O(N log N)高4D 复数数组焦栈批量生成、实时预览相位斜坡叠加每视角乘线性相位 累加中中等价于频域移动便于向量化选型建议课程设计或算法验证优先时域平移加和代码短、可解释性强要做批量焦栈或嵌入式部署优先频域方案。还有一种更工程化的做法先对光场做 4D FFT按不同 alpha 重复切片这部分在 Matlab 里用数组索引和 interpn 就能搞定后面第 4 章会给可运行版本。3.4 重聚焦系数 b 与景深换算不要对着公式硬调很多人在 alpha 和 b 之间反复试错原因是没有建立像素位移到深度之间的标定关系。常见做法是光场标定得到基线长度 B米和微透镜焦距 f米那么深度 d 对应的子孔径位移系数b B * f * (1/d_ref - 1/d)在不确定标定参数时不直接输深度输 b 的浮点范围 [-2, 2]观察焦点前后变化。演示阶段用 b 更直观b0 是参考焦点b 绝对值越大焦点偏离参考面越远。深度单位只有等标定完成后再套公式换算。4. 两种重聚焦算法的 Matlab 实现与参数选择4.1 时域平移加和用 imtranslate 写一个可直接复用的函数下面的函数直接接受子孔径图像堆和一个重聚焦系数 b输出重聚焦后的彩色图。它做了三件在真实项目里必须处理的事单精度存储省内存、边界填空值、视角数归一化。function focused ShiftAndAddRefocus(lf, b) % ShiftAndAddRefocus: 时域平移加和重聚焦 % lf: [nT, nS, nV, nU, 3] 子孔径图像堆 % b: 重聚焦系数0 表示参考焦点正负分别对应后景/前景 % 输出 focused: [nT, nS, 3] 单精度 [nT, nS, nV, nU, ~] size(lf); lf single(lf); cx (nU 1) / 2; cy (nV 1) / 2; frames cell(nV, 1); for v 1:nV tmp zeros(nT, nS, 3, single); for u 1:nU img squeeze(lf(:, :, v, u, :)); shiftX b * (u - cx); % x 方向像素位移 shiftY b * (v - cy); % y 方向像素位移 tmp tmp single(imtranslate(img, [shiftX, shiftY], ... FillValues, 0, OutputView, same)); end frames{v} tmp; end focused zeros(nT, nS, 3, single); for v 1:nV focused focused frames{v}; end focused focused / (nU * nV); end代码逻辑分三层内层循环逐视角平移累加中层把同一行视角的中间结果存到 cell外层做最终归一化。把 frames 存下来而不是直接在 parfor 里累加是为了避免 Matlab 并行池在 reduction 变量上的类型不一致问题。参数说明shiftX b * (u - cx)中 u-cx 是视角相对中心的整数偏移b 控制每单位视角偏移对应的像素数。b1 意味着相邻视角间移动 1 个像素。OutputView, same保持输出尺寸和输入一致如果要生成用于拼接的全景重聚焦可以改成full但边界黑边会更明显。FillValues决定移出边界的像素填充值。0 会产生黑边对后续评价函数有干扰边缘信息重要时改用replicate或镜像填充但速度会变慢。如果装了 Parallel Computing Toolbox把外层for v 1:nV改成parfor v 1:nV在 8 核机器上大约能快 3~5 倍。注意切片squeeze(lf(:,:,v,u,:))在 parfor 里会产生广播变量内存紧张时建议先把 lf 按 nV 维度拆成 cell 传进去。4.2 傅里叶切片法用 4D FFT 一次生成多张焦栈频域方法一次 FFT 后可以对任意 b 做切片复用因此批量重聚焦时优势明显。这里实现一个工程上更稳的相位斜坡版本对每个视角图在频域乘线性相位等价于时域亚像素平移避免了 imtranslate 的插值误差。function focusStack FFTRefocusStack(lf, bList) % FFTRefocusStack: 傅里叶域重聚焦一次生成多个 b 的重聚焦图 % lf: [nT, nS, nV, nU] 灰度光场彩色需逐通道处理 % bList: 重聚焦系数数组如 -2:0.2:2 % 输出 focusStack: [nT, nS, numel(bList)] [nT, nS, nV, nU] size(lf); % 空间频域坐标 fy (0:nT-1) - nT/2; fx (0:nS-1) - nS/2; [FX, FY] meshgrid(fx, fy); % 角度坐标相对中心 au (0:nU-1) - (nU-1)/2; av (0:nV-1) - (nV-1)/2; % 对角度维做 FFT相当于把空间维的每个像素变换到频域表达 F fft(fft(lf, nU, 4), nV, 3); focusStack zeros(nT, nS, numel(bList)); for k 1:numel(bList) b bList(k); acc zeros(nT, nS); for v 1:nV for u 1:nU % 线性相位exp(-2πi * b * (au*FX av*FY) / N) phase exp(-2i * pi * (b * (au(u)*FX av(v)*FY)) / nS); acc acc F(:, :, v, u) .* phase; end end focusStack(:, :, k) abs(acc) / (nU * nV); end end这段代码先把光场在角度维做 FFT再在空间频域对每个视角乘相位斜坡。原理上时域的整数像素平移对应频域的线性相位所以 b 不再要求是整数亚像素精度天然支持。相比直接做 4D 切片这个写法内存占用更低循环结构也更清晰。参数说明bList控制输出焦栈的深度范围和步长。步长太小导致相邻两帧差异不明显步长太大会跳过焦点平面。我一般先按 0.2 步长扫一遍找焦点大致范围再在焦平面附近用 0.05 步长细化。相位公式里除以 nS 是让空间频率归一化保证 b 的物理含义和时域版本一致单位为像素/视角。该实现只针对灰度图彩色图需要对三个通道分别调用或把lf扩展成 5 维后对第 5 维循环。4.3 参数怎么定视角数、b 的范围、内存预估视角数 nU×nV 是光场最重要的参数。仿真数据里 5×5 足够验证算法真实 Lytro 数据一般有 14×14 以上视角但边缘视角光圈渐晕严重实际使用时常裁剪成 7×7 或 9×9。视角越多重聚焦后分辨率越高但角度采样率越高对内存和计算量的压力也越大。视角数光场内存512×512×3单帧重聚焦耗时时域适用场景5×5约 37 MB约 80 ms算法验证、课设演示7×7约 73 MB约 150 ms平衡精度与速度9×9约 121 MB约 250 ms生成高质量焦栈上表按 single 精度估算实际中如果光场存成 double内存翻倍8 GB 内存的机器跑 9×9 会很紧张。建议在生成光场或读入数据后立刻single()转换这个习惯比任何“优化工具箱”的设置都直接有效。b 的范围取决于场景视差幅度。仿真数据里视差范围 [0,2] 对应像素位移 0~40b 取 [-2,2] 基本覆盖真实场景要先计算最大视差取一张子孔径图上的明显特征点在另一视角子孔径图里找对应点像素距离就是该深度的最大位移b 的范围设为这个最大位移的 ±1.5 倍即可。4.4 真实数据验证从解码到重聚焦的完整流程拿到 Lytro 或 Raytrix 的真实数据时解码部分各厂商差异很大。这里只给通用的后续流程骨架具体解码函数需要按相机 SDK 替换白平衡用白图像对微透镜图做逐像素增益校正消除微透镜渐晕和传感器暗角。微透镜中心标定对白图像做阈值分割拟合每个微透镜中心的整数坐标。按第 2.2 节重构子孔径图像堆。对子孔径图做亮度归一化每张子孔径图除以自身均值消除角度维亮度不均。调用第 4.1 或 4.2 节的重聚焦函数先用大步长 b 找到焦点。真实数据和仿真数据的最大区别是存在渐晕、微透镜畸变和颜色串扰因此在重构后最好先做一次“中心视角与边缘视角的对齐检查”把两个视角的图叠加显示成红青图观察物体边缘是否有颜色拖影有拖影说明旋转角或者 pitch 没标定准。5. 用 App Designer 给重聚焦算法加一个可拖动的演示界面5.1 最小可运行界面一个滑块、一张图、一次回调演示阶段最忌讳每次调参都回命令行改数字。我一般用 App Designer 的 uifigure 搭配 uislider 和 uiaxes 做一个极简界面总代码量不到 40 行。下面这个脚本可以直接复制运行前提是工作区里已有 lf 变量。function RefocusDemo(lf) % RefocusDemo: 拖动滑块实时重聚焦演示 % lf: [nT, nS, nV, nU, 3] 子孔径图像堆 fig uifigure(Name, 光场相机重聚焦演示, Position, [100 100 800 700]); g uigridlayout(fig, [2 1]); g.RowHeight {1x, 60}; ax uiaxes(g); slider uislider(g, Limits, [-2 2], Value, 0); slider.MajorTicks -2:0.5:2; slider.ValueChangedFcn (s, ~) updateFocus(s.Value); cache single(lf); % 回调里不再读原始变量避免类型转换重复执行 isBusy false; function updateFocus(b) if isBusy return; end isBusy true; drawnow limitrate; img ShiftAndAddRefocus(cache, b); imshow(ax, im2uint8(img)); title(ax, sprintf(重聚焦系数 b %.2f, b)); isBusy false; end end回调逻辑很直接滑块数值变化时调用 updateFocus内部用 isBusy 标志做防抖。拖动滑块时 ValueChangedFcn 触发频率远高于屏幕刷新率如果不防抖回调会堆积导致界面卡死drawnow limitrate再配合标志位让界面保持流畅。运行前确认ShiftAndAddRefocus函数在 Matlab 路径下。如果要调 b 的实时响应速度优先缩减 lf 的空间分辨率而不是视角数——空间分辨率对视觉影响不大但能把单帧重聚焦时间降一半以上。5.2 多焦点对比视图近景、中景、远景三张图同屏输出只给一个 b 的实时画面观众很难感受到“重聚焦”和“缩放”的区别。我习惯在界面下方加三个固定 b 值的缩略图分别取前景、参考面、背景一眼就能看出不同深度的清晰度变化。% 在 uigridlayout 中增加一行放三个 uiaxes bg uigridlayout(g, [1 3]); previews gobjects(1, 3); bList [-1.5, 0, 1.5]; labels {前景聚焦, 参考面, 背景聚焦}; for k 1:3 previews(k) uiaxes(bg); img ShiftAndAddRefocus(cache, bList(k)); imshow(previews(k), im2uint8(img)); title(previews(k), labels{k}, FontSize, 6); end如果希望对比时能看到清晰度差异的细节可以在每张图上叠加一个插值放大的局部窗口用xlim/ylim锁定同一 ROI 区域。这个技巧在答辩演示时很有用因为投影仪上 512×512 的图很难看出微小差异局部放大才是判断焦点是否准确的依据。5.3 演示完怎么导出结果图像和矢量图两种出口界面里的图可能还需要写进报告或论文。uifigure 和传统 figure 的导出方式不同不要用旧版的print(gcf, -dpng)对 uiaxes 截图涉及 4K 屏时容易出现空白。推荐用exportgraphics导出分辨率由第二个参数控制。% 导出当前重聚焦图像 exportgraphics(ax, refocus_result.png, Resolution, 300); % 导出论文矢量图EPS fig2 figure(Visible, off); imshow(im2uint8(img)); exportgraphics(fig2, refocus_result.eps, ContentType, vector); close(fig2);导出 EPS 时先在不可见 Figure 里重新绘制一遍避免界面组件混入图像。想进一步压缩 EPS 文件大小可以先转 SVG 再用 Inkscape 转 EPSMatlab 直接输出的 EPS 体积通常偏大。exportgraphics从 R2020a 起可用旧版本只好退回print(gcf, -depsc2)但透明背景和中文标题的支持要看具体环境。6. 聚焦质量怎么验证评价函数与 3 个高频坑6.1 用自动评价函数找最佳 b而不是靠肉眼拖动滑块找到的 b 值可以靠眼睛确认但要定量描述或多帧对比时得用一个与清晰度正相关的数字把这个数在对焦过程中最大化。图像领域称这一类为聚焦评价函数同一个评价函数还可以在光场深度估计里直接当置信度用。下面是用 Matlab 的 image processing 能力实现三种常用评价函数。function score FocusMeasure(img, method) % FocusMeasure: 单帧图像清晰度 % img: uint8 或 double彩色图自动转灰度 % method: variance / lap / tenengrad g double(rgb2gray(img)); switch method case variance score var(g(:)); case lap h fspecial(laplacian, 0.2); score mean2(imfilter(g, h, replicate) .^ 2); case tenengrad sx fspecial(sobel); gy imfilter(g, sx, replicate); gx imfilter(g, sx, replicate); score mean2(gx .^ 2 gy .^ 2); end end配合一个循环就能画出聚焦评价曲线峰值所在的 b 就是最优焦点bRange -2:0.05:2; scores zeros(size(bRange)); for k 1:numel(bRange) img ShiftAndAddRefocus(cache, bRange(k)); scores(k) FocusMeasure(img, lap); end [~, idx] max(scores); fprintf(最优 b %.2f\n, bRange(idx));三种评价函数对噪声的敏感度不同灰度方差对平坦区域很钝拉普拉斯能量对边缘锐度敏感但怕噪声Tenengrad 用 Sobel 梯度所以最稳。实际使用中推荐拉普拉斯或 Tenengrad如果重聚焦图边界有黑边干扰裁剪中心 80% 区域再算否则评价结果会被黑边主导。6.2 三个高频坑中心偏移、边缘黑边、频域混叠重聚焦在真实数据上跑不干净绝大多数是这三个原因之一。微透镜中心偏移。重构子孔径图时如果微透镜中心坐标偏了半个像素叠加后的图像会叠加一层虚影。判断方法用中心视角和相邻视角做差观察残差图里是否有一侧亮的边缘拖影。修复方案是在白图像上做亚像素质心拟合不要直接用整数网格坐标。边缘黑边。平移加和时FillValues0会在图边缘产生黑边影响评价函数。常见做法是重聚焦后裁剪到有效区域平移量最大为 |b|×(nU-1)/2所以有效区域四边各裁剪这个像素数。另一种思路是对光场先做镜像扩边再重聚焦完事后再裁回原尺寸边缘质量会好很多。频域混叠。傅里叶切片法在视角数不足时重聚焦图上会出现周期性的波纹或重影尤其在高频纹理区域明显。缓解手段有三个重聚焦前对光场做低通滤波对角度维补零到 2 倍后再做 FFT或者在相位斜坡叠加时对 b 做亚像素步进比如 0.05混叠会随步长细化而减轻但不能完全消除。时域平移加和在同样条件下不会有这类混叠这是它在精度敏感场景下仍然被选用的原因。验证整个流程是否正确的最后一步是把 b 从 -2 到 2 逐帧播放观察焦点是否从前景平滑移动到背景且评价曲线的峰值位置与视觉清晰面重合。如果峰值位置和直观判断不一致优先检查是否有边缘黑边干扰评价函数其次检查光场的角度维顺序是否反了——lf(:, :, v, u, :)和lf(:, :, u, v, :)在调用ShiftAndAddRefocus时是不同结果这种维度顺序错位在代码审查里很难发现建议在函数入口加一行size(lf)打印确认。本文还有配套的精品资源点击获取