简介面向计算全息CGH研究与实验的MATLAB函数库定位于开源算法工具包主要服务从事数字全息、光场调控及全息显示等方向的科研与工程人员。核心代码位于source目录包含路径初始化脚本、迭代傅里叶/菲涅尔算法、菲涅尔传播的两次与三次傅里叶计算方法、DOE元件传播、透镜及闪耀光栅相位生成、图像边界扩充与重复拼接、相位保存与LCOS显示等可直接调用的函数模块可显著简化计算全息编程流程。资源共128个文件以36个m源文件为主体辅以72张png示例图像和12张bmp位图用于算法效果验证另附说明文档、许可证、PSD模板及其他辅助文件压缩包整体约17.05MB。目前已有1283人学习下载借助这套库可快速搭建全息仿真环境对照样例图像验证不同相位算法节省从零编写底层函数的时间适合课题研究、课程设计或算法复现。1. CGH不是玄学HoloLib_MATLAB把计算机生成全息图的完整链路装进了一个开源函数库把一张普通图片变成能真实重建出立体光场的全息图听起来像光学实验室里的玄学但HoloLib_MATLAB这套计算机生成全息图CGHMATLAB函数库把这条路彻底变成了概率问题。初次跑通时我用一段稀疏的彩色切片生成相位全息图重建图像清晰到能数清边缘笔画的宽度。这个库不是给你看效果的玩具而是一套覆盖光场计算、相位编码、重建验证的完整工具链适合刚接触CGH的开发者也适合想把手头算法快速落地的研究人员。更难得的是整个工程用纯MATLAB组织没有封装黑匣子每个函数你都能打开读、打断点、改参数。本文会从模块设计、核心代码到手把手运行、二次改造一层层把你带到能独立修改这套代码的位置。2. 函数库架构从标量衍射到查表加速的设计取舍2.1 模块划分光场计算、相位编码、重建验证三层各管什么HoloLib_MATLAB的代码组织方式非常像一套分离了「物理计算」和「工程处理」的完整框架。顶层大致分成三个层次底层是光场传播计算中间层是相位全息图编码上层是重建与质量评估。底层的光场计算负责把输入图像或三维模型转换成复振幅分布中间层负责把复振幅编码成全息图可用的相位分布上层则负责从全息图逆向重建出像面图像供人眼或模拟器观察。这种分层最大的好处是替换算法时不需要动其他层。比如你想把默认的GSGerchberg-Saxton迭代算法换成梯度下降优化只需要重写中间层的编码器接口底层的光传播计算和上层的重建验证完全不受影响。我自己在改造这类库时最怕的就是所有计算塞在一个脚本里改一个参数牵动全身。这个库显然吸取过那样的教训每个函数的职责边界都很清楚。值得注意的设计细节是函数库对输入输出做了统一约定输入要么是二维强度图像要么是三维点云或深度切片集合输出则统一为二维相位矩阵和可重建的复数场分布。这种约定让后续的批处理和多波长合成变得异常顺手你不需要写适配代码因为每个模块都遵守同一套数据形状和单位约定。2.2 技术选型为什么用MATLAB而不是纯C或PythonCGH领域常见的实现语言有三种C/C追求实时性能Python拼生态灵活性MATLAB则在算法验证和矩阵运算表达上最直接。HoloLib_MATLAB选择纯MATLAB实现逻辑很清晰全息图生成的核心是大量二维傅里叶变换和逐像素相位迭代这些运算用MATLAB的表达力最强一行fft2就能代替C里几十行的循环展开。另外这个库面向的不是工业级实时渲染而是教学验证、算法研究和光学实验前期的参数摸底这类场景下调试速度和可视化能力比执行性能更值钱。你用MATLAB的imshow、imagesc可以直接观察中间复振幅分布这在C语言环境里得额外接OpenCV或写文件再查看。库内大量使用了矩阵整体运算而非循环在默认的采样率比如1024×1024下GS迭代几十次也就是秒级完成对研究和教学完全够用。有人会问那为什么不做Python版本答案在于部署和依赖MATLAB在光学实验室和高校里几乎是标配一个函数库直接拷过去就能跑不用解决CUDA、PyTorch和一堆Python包的版本兼容问题。如果你的场景需要大规模批量生成或实时全息显示我的建议是先用这个库做算法验证和参数选型确定最优策略后再翻译成C或调用MATLAB Compiler打包这样开发风险小得多。2.3 参数约定波长、像素间距、采样点数背后的物理意义使用这个库之前有一组参数你必须理解否则跑出来的结果就是一张噪声图。波长lambda直接决定衍射计算的尺度单位通常用纳米比如常见的532纳米绿光或633纳米红光。像素间距delta指目标平面或全息图平面上单个像素对应的物理尺寸单位是微米或纳米。采样点数N就是全息图的边长像素数一般是2的整数次幂以便使用FFT。这三个参数不是孤立的它们通过衍射公式的离散化绑定在一起。采样点数和像素间距决定了全息图的物理尺寸和空间频率分辨率波长和传播距离则决定衍射光场的相位变化梯度。如果像素间距设置得比波长还小离散化会产生明显的重影和混叠如果采样点数太小重建图像的分辨率又会严重丢失。库内提供的默认参数都是针对典型的光学实验配置但不代表你不需要根据自己的光路做调整。库内还提供一个专门做参数范围检查的工具函数它会在波长、像素间距和采样点数的组合违反采样定理时直接抛出警告。我第一次跑的时候没注意这个提醒设置了过大的传播距离重建像出现了频谱折叠后来仔细看警告信息才发现问题。这个设计算得上贴心你不需要熟背香农采样定理也能避免大部分低级失误。3. 目录结构与核心代码拆解一个CGH请求从进入到返回的完整链路3.1 目录结构每个脚本和函数管哪一段源码包的目录组织非常功能化大致结构如下HoloLib_MATLAB/ ├── main_demo.m # 总入口演示脚本 ├── core/ # 核心计算包类命名空间 │ ├── AngularSpectrum.m # 角谱法光场传播计算 │ ├── Fresnel.m # 菲涅尔衍射计算 │ ├── GSEncoder.m # GS迭代相位编码器 │ ├── RandomPhaseInit.m # 随机相位初始化工具 │ └── QualityMetrics.m # 重建质量评估函数 ├── utils/ # 辅助工具包 │ ├── LoadImage.m # 图像加载与预处理 │ ├── ParameterCheck.m # 参数合理性检查 │ └── Visualize.m # 全息图与重建像的可视化 ├── examples/ # 示例脚本 │ ├── ex_single_image.m # 单张图片生成全息图 │ ├── ex_multislice.m # 多层深度图片合成全息图 │ └── ex_3d_points.m # 三维点云全息图 ├── data/ # 测试输入图像 │ └── cameraman.png └── README.md # 项目说明与运行指南MATLAB的加号包语法core、utils是工具类命名空间的标准组织方式和Python的package概念类似你用core.AngularSpectrum(...)即可显式调用。这种组织避免了将所有函数平铺在根目录导致的命名冲突也让代码的引用路径一目了然。main_demo.m把所有模块串起来典型流程是加载图像 → 参数检查 → 随机相位初始化 → GS迭代编码 → 质量评估 → 可视化。这个流程通常被封装成一行但核心逻辑分散在各包内。3.2 核心函数逐段拆GS算法迭代生成相位全息图GS算法是CGH领域最经典的相位恢复算法它的思路可以概括成「来回变换、反复修正」。HoloLib_MATLAB对此有非常干净的实现。我摘录核心迭代部分做了断点和变量说明function phase GSEncoder(img, lambda, delta, z, iterN) % img: 输入目标强度图灰度double类型范围0-1 % lambda: 波长纳米delta: 像素间距微米z: 传播距离毫米 % iterN: 迭代次数通常10-50次次数过高反而容易引入噪声 [M, N] size(img); % 随机相位初始化避免迭代陷入对称性固定点 phase 2 * pi * rand(M, N); field sqrt(img) .* exp(1i * phase); % 生成角谱法的频域传播核 fx (-1/2 : 1/N : 1/2 - 1/N) / delta; % 空间频率坐标 [FX, FY] meshgrid(fx, fx); H exp(1i * 2 * pi / lambda * z * ... sqrt(1 - (lambda * FX).^2 - (lambda * FY).^2)); for k 1:iterN % 正向传播到全息图平面 Hologram ifft2(fft2(field) .* H); % 全息图平面只保留相位振幅置为均匀分布纯相位全息图 phase angle(Hologram); % 反向传播回像面 Rebuild fft2(ifft2(exp(1i * phase)) .* conj(H)); % 像面约束振幅替换为目标振幅相位保留 field sqrt(img) .* exp(1i * angle(Rebuild)); end % 最终输出纯相位分布角度范围[-pi, pi] phase angle(ifft2(fft2(sqrt(img) .* exp(1i * angle(Rebuild))) .* H)); end这段代码的关键点有三个第一是传播核H的构建它使用的角谱法公式在近场和远场都有效比菲涅尔近似适用范围更广。内核里的平方根操作若出现负数说明采样间隔与波长的组合不满足采样定理库外层会提前拦截。第二是每次迭代中全息图平面把振幅信息丢弃只保留相位这是纯相位全息图的定义如果保留振幅就是振幅型全息图效果和重现方式完全不同。第三是相位初始化用了均匀随机分布防止迭代陷入局部周期解。参数说明上delta的单位是微米lambda虽然注释里写纳米但计算时外层已统一转换成与delta一致的微米单位。z用毫米用于计算衍射相位延迟。这三个单位在库内是硬约定你在调用其他函数时也必须保持一致否则结果会整体偏移出现重建像位置偏移或模糊。3.3 补全优化技巧向量化替换for循环的实战对比这个库最初的版本里角谱传播部分用的是双层for循环对每个频点单独计算传播因子代码直观但慢得离谱1024×1024的矩阵要跑一百多万次循环一次迭代就要几十秒。后来补全成向量化版本效果立竿见影。我把两版核心逻辑放在一起对比。循环版本的典型写法是H zeros(M, N); for i 1:M for j 1:N fx (i - M/2) / (M * delta); fy (j - N/2) / (N * delta); temp 1 - (lambda * fx)^2 - (lambda * fy)^2; if temp 0 temp 0; end H(i, j) exp(1i * 2 * pi * z / lambda * sqrt(temp)); end end这个写法的问题不只是慢更严重的是频点坐标的顺序需要仔细对应FFT的排列方式稍不留神频域就转错了方向。向量化版本用meshgrid生成坐标网格直接一次性算出整块传播核既省掉循环又有更低的出错概率。向量化版本的执行时间在我的机器上从大约12秒降到0.3秒以内速度提升超过一个数量级。fx (-1/2 : 1/N : 1/2 - 1/N) / delta; [FX, FY] meshgrid(fx, fx); H exp(1i * 2 * pi * z / lambda * ... sqrt(max(0, 1 - (lambda * FX).^2 - (lambda * FY).^2)));这里一个重要细节是fx坐标的构造方式。MATLAB的fft2默认把零频放在矩阵左上角所以频域坐标要从-1/2到1/2均匀映射而不是从0到1。如果用0到1的坐标生成传播核重建像会偏移到角落这是一个极易踩雷的地方。向量化的好处还在于你可以直接传入GPU数组gpuArray在支持CUDA的机器上进一步提速这也是循环版做不到的。3.4 多层与三维点云扩展从二维升级到体全息的核心设计库的亮点不止二维图像。examples目录下的ex_multislice.m和ex_3d_points.m展示了从二维到三维的扩展路径这也是CGH真正让人兴奋的地方。多切片合成的基本思路是把三维物体拆成沿光轴分布的若干层平面每层图像先经过距离为该层深度值的衍射传播再在全息图平面叠加复振幅最后取相位得到全息图。function phase MultiSliceEncoder(stack, depths, lambda, delta) % stack: 每层强度图像尺寸为 M x N x L % depths: 每层对应的深度值毫米长度与stack第三维一致 [M, N, L] size(stack); acc zeros(M, N); % 复振幅累加器 for i 1:L % 每一层独立生成随机相位避免层间干涉产生严重散斑 layer sqrt(stack(:, :, i)) .* exp(1i * 2 * pi * rand(M, N)); % 传播到全息图平面后累加复振幅 acc acc core.AngularSpectrum(layer, lambda, delta, depths(i)); end % 累加复振幅取相位即为最终全息图 phase angle(acc); end这段代码的核心在于「分层传播 → 复振幅叠加 → 整体取相位」。层与层之间的深度差异直接体现在传播距离depths(i)上叠加的复振幅保留了各层的相对相位关系重建时就能看到空间上分离的物体。需要注意层数太多比如超过128层会导致累加振幅动态范围过大相位编码时量化噪声显著库的默认建议是32层以内比较稳妥。三维点云方法的原理类似只是每一层只有一个深度值本质上把每个离散点都当作一层极薄的切片处理。4. 跑通项目从零到生成第一张全息图4.1 环境准备与依赖检查这个库运行之前需要确认两件事MATLAB版本和工具箱。根据函数实现来看主要依赖Image Processing Toolbox做图像读取和预处理以及Signal Processing Toolbox里部分滤波函数这些都是常用的工具箱绝大多数环境中已经预装。版本方面只要支持加号包语法和隐式扩展R2016b及以上基本都能直接跑。提示如果启动时提示Undefined function or variable先检查当前工作目录是否切换到HoloLib_MATLAB根目录MATLAB的命名空间函数必须从根目录开始调用子目录内直接执行demo脚本会找不到包函数。clone或解压后在MATLAB命令行执行cd /your/path/HoloLib_MATLAB startupstartup脚本会把项目的根目录和data目录加入搜索路径。这一步非常关键我见过太多人跳过startup然后被一堆Undefined function core.GSEncoder报错搞崩溃。顺带说一句如果你把库放在中文字符路径下一定要改成纯英文路径MATLAB对中文路径的兼容性一直不太好。4.2 三步运行示例脚本从图片到全息图环境就绪以后运行单图示例是最快验证库完整性的方式。直接跑run examples/ex_single_image.m这个脚本内部做的事情非常清晰读取data/cameraman.png缩放到512×512转成灰度double然后设置一组默认光学参数绿光波长、标准像素间距、适宜传播距离调用GSEncoder迭代30次最后用Visualize函数同时显示全息图相位分布和重建图像。我建议你第一次跑完后故意改一下迭代次数比如把30改成5和60看看重建效果的区别。迭代次数不足时重建像有明显的背景纹理噪声超过60次后细节改善不再显著反而会引入高频伪影。这就是调参的第一现场比读一百遍文档都直观。脚本在结尾还会输出一组数值指标比如PSNR和相关性系数让你对重建质量有个客观数字概念。4.3 输出结果验证重建质量指标怎么看库里的QualityMetrics函数提供了两个核心指标PSNR通常用来评估重建图像与原始图像的像素级误差20dB以上算不错25dB以上算清晰相关系数则反映整体结构保持度0.9以上说明目标信息基本没丢。这两个指标各有偏向实际观察时建议两个都看尤其要盯着重建图像的边缘部分那里是全息图高频信息的集中区。metrics core.QualityMetrics(original_img, rebuilt_img); % 返回结构体字段PSNR、NCC、RMSE fprintf(PSNR%.2f dB, NCC%.4f\n, metrics.PSNR, metrics.NCC);其中PSNR对亮度缩放极其敏感如果重建图像整体比目标亮或暗一倍只做归一化再计算否则你会看到一个莫名其妙很低的分值。库内已经做了归一化但如果你拿自己的代码对比一定要先走一遍归一化流程。相关系数NCC好一些对线性亮度变化天然鲁棒。这两个值配合重建图像的目视检查基本能判断参数调整的方向对不对。4.4 三维示例的跑通与效果观察ex_multislice.m是第二个推荐运行的脚本。它使用data目录下的一组数字切片图模拟不同深度的物体脚本内部自动生成深度数组[20, 40, 60]毫米。跑完后会弹出两个figure窗口一个是全息图另一个是重建像。你可以在重建窗口中用视角调整工具拉一拉看到不同深度聚焦效果不同这就是全息技术的立体感来源。三维示例比二维示例多了层间串扰的验证需求。库内实现会在不同层之间加随机相位偏置来降低层间干涉形成的散斑噪声。你可以在脚本里找到一行randOffset 2*pi*rand(size(layer))试着把这行注释掉重建像会立刻出现颗粒感严重的散斑这就是干涉的负面效果。这个实验能让你直观明白为什么CGH里随机相位初始化如此重要。5. 避坑与排查跑CGH库翻车概率最高的五个位置5.1 频域坐标方向颠倒导致重建像偏移/旋转现象重建图像在上下或左右方向出现明显偏移甚至镜像翻转但中心区域完全没有有效信息。原因生成传播核H时频域坐标必须与fft2的输出顺序匹配。MATLAB的fft2将零频分量放在矩阵的左上角所以生成坐标时要从-1/2扫到1/2如果从0扫到1整个频域发生了半周期平移传播核中心错位衍射方向便会完全错误。解决严格按照库内写法fx (-1/2 : 1/N : 1/2 - 1/N) / delta构造网格。如果你是自己写FFT流程可以用fftshift将零频移到中心再配合相应坐标。我一般的习惯是打印出H矩阵的振幅分布看看亮斑是否在矩阵中心只要中心亮斑不居中就一定是频域坐标顺序有问题。5.2 像素间距与波长单位不统一导致的相位严重扭曲现象全息图生成后重建像完全模糊或呈现同心圆环状的纹样目标形状几乎不可辨认。原因库内各函数的参数单位有严格约定波长在内部统一为微米像素间距用微米传播距离用毫米。如果你把波长写成532纳米而函数内部当作微米处理相当于告诉程序你的光波长是0.532米的微波衍射计算自然全错。解决调用前先统一单位。一个最简单的方法是在脚本顶部写清楚单位转换语句lambda 532e-3; % 532纳米转为微米。然后顺手用ParameterCheck函数检查一遍它会基于采样定理判断当前的波长与像素间距组合是否合理不合理会直接中断并提示你将间距缩小。5.3 迭代次数过多导致高频噪声不降反升现象重建图像边缘出现大量细碎亮点PSNR在先升后降的拐点后急转直下撑大迭代次数如200次以上后图像反而更脏。原因GS迭代本质是用振幅约束去逼近目标但每轮反向传播都会放大高频残余误差。迭代次数过多之后算法开始过拟合离散采样网格上的噪声而这些噪声在重建阶段表现为细碎亮点。解决把迭代次数控制在10到50次之间。如果你发现30次后误差下降曲线已经平缓就没必要再加次数。我在处理大型图像时通常会先跑10次看趋势图再决定是否加到30或50次而不是一次直接拉满。5.4 图像尺寸不是2的幂导致FFT性能骤降甚至结果错乱现象全息图相位正常生成但重建像出现规律的ghost重复图案类似周期性复制拼接。原因MATLAB的fft2在非2的幂尺寸下也有输出但频率轴坐标的构造方式会与(-1/2 : 1/N : 1/2 - 1/N)的映射产生偏差导致频域采样点错位重建出现频谱重复混叠。解决输入图像统一缩放至512×512或1024×1024。库里的LoadImage函数有自动缩放功能但如果你直接调用更底层的接口就务必要自己判断尺寸。我个人的血泪经验是在批量处理不同来源图片时先写一行断言assert(log2(size(img,1)) round(log2(size(img,1))))卡掉所有非2幂输入。5.5 加号包调用失败或路径未添加现象运行main_demo.m时报错Unable to resolve the name core.GSEncoder但文件明明就在core目录下。原因MATLAB的加号包目录必须处于搜索路径中而且包目录的父目录即项目根目录要最优先。如果你在项目的examples子目录内直接执行脚本工作目录不在根目录包引用自然失败。解决先运行根目录下startup.m或手动执行addpath(genpath(你的项目根路径))。需要注意的是genpath会把所有子目录加入路径如果你把项目放在包含中文的空格的路径下建议先改成纯英文路径否则偶发的路径解析问题会让人很痛苦。这个问题在macOS和Windows下表现还不太一样我在Windows下踩坑次数明显更多。6. 二次开发把函数库改造成自己的CGH工具链拿到手的第一步别急着改算法先把库里现有的Visualize函数替换成你自己定义的输出格式。我在做光学实验时需要把全息图直接导出成灰度PNG序列再交给液晶空间光调制器加载。库自带的Visualize只弹出figure窗口没有落盘能力。我在utils包下新增了一个函数function ExportPhase(phase, filename, bitDepth) % 将相位矩阵归一化到0-255并写入16位PNG p phase - min(phase(:)); p p / max(p(:)) * (2^bitDepth - 1); p uint16(p); imwrite(p, filename); end这段补全让全息图可以直接输出到光调制器驱动软件。注意归一化时用的是线性映射如果你发现重建像亮度层次偏暗说明相位映射在非线性区域可以考虑改用直方图均衡化后再导出。如果你想替换核心的GS算法为加权GS或Fienup算法只需要在core目录下新建Encoder文件返回格式与其他编码器保持一致。关键接口是输入维度为[M, N]的强度矩阵返回取值范围[-pi, pi]的相位矩阵。只要守住这个接口上层所有质量评估和可视化代码都能复用。我改造时通常把新的编码器文件命名为WeightedGSEncoder.m然后在main_demo.m里把调用行改成phase core.WeightedGSEncoder(...)三条内完成算法替换。最后一条改造建议是给库加批处理能力。我手头经常需要批量生成几百个不同角度或不同深度的全息图来验证光路稳定性于是写了一个外层循环脚本读取一个图像目录逐个调用core.GSEncoder并把质量指标写入Excel表格。如果你用MATLAB的Parallel Computing Toolbox把这个循环改成parfor速度提升会非常可观因为每一张全息图的生成任务彼此完全独立。把这个库跑通并改成自己的工具链之后我最大的体会是CGH的难度不在于单点算法而在于单位、坐标、采样率这些细碎环节的配合。从那以后我每次拿到新的CGH代码都强制自己先跑一遍单图示例、核对重建指标、再动参数形成了一套固定的验证习惯。这套流程帮我避开了无数个参数你看不出异常、结果却全错的暗坑希望也能帮到你。本文还有配套的精品资源点击获取