企业官网建设流程全解析

简介：一套基于MATLAB的数字全息仿真代码包，面向光学成像、计算机图像处理方向的研究者和学生，可用于计算全息图的制作与再现实验。包内共18个文件，以6个.m脚本为主体，覆盖菲涅尔衍射全息图生成、在线全息重建、并行处理等典型算法；同时附有2个jpg、2个bmp、1个tif、1个png共6张示例图像，包括经典cameraman测试图，以及docx、wps说明文档，便于对照代码理解每一环节的输出效果。算法代码从全息记录面的干涉场计算，到利用傅里叶变换进行波前重建，形成了完整的数字全息处理流程；资源中还包含多个变体和备份脚本，可观察不同参数设置对再现图像质量的影响，适合在修改衍射距离、波长等条件后进一步分析。整个压缩包约515KB，体量虽小但信息完整，已有690人学习下载，通过这份资料不仅能跑通全息图生成到再现的完整示例，还能借鉴代码结构设计自己的仿真实验，为深入全息三维显示、防伪识别等应用打下基础，是入门数字全息编程和验证再现算法的实用参考。

1. 数字全息再现：zip 文件包里藏的科学，比压缩算法更值钱

数字全息.zip_cryni1_全息_全息再现_全息图再现这个文件名乍看像随手存的压缩包，实际是一个含噪的采集结果：cryni1像是采集批次或样品编号，全息出现了三次，说明这批数据是同一个样品在不同再现条件下的反复处理。数字全息要解决的核心问题不是存储格式，而是从一张灰度干涉图（全息图）里恢复出物光波的振幅和相位，进而算出微米甚至纳米级的三维形貌。这个过程的专业说法叫全息再现，最常用的算法是菲涅尔衍射积分和角谱法。读这篇文章的人，要么手头正握着类似命名的 .zip，要么在用显微镜做相衬观察却不知道如何从干涉条纹里挖出定量数据，要么在排查为什么"明明拍了图却什么都没有重建出来"。下面从头到尾讲一遍数字全息图的读取、预处理、数值衍射复现和相位提取，代码直接用 Python 写，参数全部给到能跑通为止。

2. 数字全息图的读入与预处理：把干涉条纹变成可计算的复数场

2.1 先认清手上有什么：参考光、物光和全息图的角色

任何一个数字全息记录系统，离不开三样东西：物光 O、参考光 R 和记录面上的强度分布 H。CCD 或 CMOS 传感器只能记录光强，所以干涉条纹被写成

H(x, y) = |O|² + |R|² + O·R* + O*·R

前两项是直流项（零级像的源头），第三项是物光与共轭参考光的干涉项（包含原始物光信息），第四项是孪生像。全息再现的目标，就是把第三项从 H 里分离出来，并让它在数值上重新传播回物平面。

解压数字全息.zip之后，通常会看到 R、O、H 三张图，或直接看到一张带斜条纹的灰度图。cryni1这种后缀一般对应采集时的激光波长或曝光参数标识，读数据前先确认像素位深（常见 8bit、16bit），再确认图像尺寸。位深错了，后面的强度归一化和频谱滤波全都白做。常见做法是先用 PIL 或 imageio 读入，看最小最大值，判断是否需要除以 4095（12bit）或 65535（16bit）。

2.2 减背景与归一化：直流项必须尽可能压平

直接对原始全息图做再现，DC 项会在再现像中心形成一个巨大亮斑，把细节全淹没。标准的预处理分三步：暗场减除、平场校正、整体归一化。暗场（不开激光时采集）减掉传感器的固定噪声；平场（无样品时采集参考光）用来校正照明不均匀。

import numpy as np from imageio.v2 import imread from scipy.fft import fft2, ifft2, fftshift hologram = imread('cryni1_H.tif').astype(np.float64) dark = imread('cryni1_dark.tif').astype(np.float64) flat = imread('cryni1_flat.tif').astype(np.float64) # 结构：先减暗场，再除以平场，最后做零均值归一化 H = (hologram - dark) / (flat - dark + 1e-6) H -= H.mean() H /= H.std()

这段代码里的flat - dark对应有效照明分布，1e-6是防止除零。归一化到零均值、单位标准差之后，DC 项被整体压低，后续频谱里三个峰（零级、正一级、负一级）对比度会明显提升。需要说明的是，平场校正不是必需的——如果激光照明本身非常均匀，跳过这一步也能出结果，但加了之后频谱边缘的环状伪影会少很多。

2.3 频谱滤波：从干涉图里把物光波单独剥出来

全息再现的本质是滤波：在空间频谱域把 +1 级（或 -1 级）的频谱切出来，移到中心，再做一次逆傅里叶变换，得到一个复数场。这个复数场就是传感器平面上重建出的物光波分布，幅值对应光强，相角对应光程差。

F = fftshift(fft2(H)) rows, cols = H.shape mask = np.zeros_like(F, dtype=np.float64) cy, cx = rows // 2, cols // 2 # 参数说明：peak_y、peak_x 是 +1 级频谱峰的中心位置，需从频谱强度图里目视读出 peak_y, peak_x = 178, 223 r = 36 # 滤波半径，按条纹载频的 2 倍宽度给 yy, xx = np.ogrid[:rows, :cols] mask = ((yy - peak_y) ** 2 + (xx - peak_x) ** 2) <= r ** 2 F_filtered = F * mask complex_field = ifft2(fftshift(F_filtered))

peak 位置和半径 r 是这段代码里最依赖人工判断的两个参数。半径太小会切掉高频信息，形貌边缘变钝；半径太大则把零级残留或孪生像的一部分圈进来，再现相位图上出现规律性波纹。经验做法是：先打印np.abs(fftshift(F))的对数强度图，用 matplotlib 直接看峰的位置，再用10 * np.log10(np.abs(F).max()) - 20 dB作为阈值自动找峰的范围。需要提醒的是，不要试图用一个固定半径跑完所有样品——光学平台的震动状态、样品散射强弱都会让峰的形状变化，滤波器半径跟着样品走才稳定。

提示：频谱峰不在整像素点上时，掩模边缘会产生截断伪影。可靠做法是加一个高斯渐变边沿而非硬圆形掩模，半径外 5 个像素内让权重从 1 平滑降到 0。

3. 数字全息再现的数值算法：菲涅尔卷积、角谱法怎么选、怎么算

3.1 为什么实际处理很少用直接菲涅尔变换

拿到传感器面上的复数场 u(x, y) 之后，全息再现要解决的是"让光场继续往前传"的问题。严格解是瑞利-索末菲衍射积分，计算量太大；近轴近似下退化为菲涅尔衍射积分

U(ξ, η) = exp(jkz) / (jλz) · exp[jk(ξ²+η²)/2z] · ∫∫ u(x, y) · exp[jk(x²+y²)/2z] · exp[-j2π(xξ+yη) / (λz)] dxdy

直接按这个式子离散化，输出的像素尺寸会随 z 改变，而且为了满足采样率，z 被限制在很小的范围内。这导致一个实际问题：当你对同一张全息图试 z=50mm 和 z=55mm 时，两幅再现像的横向分辨率都变了，没法直接对比。更难受的是，当 z 太小（比如显微配置下只有几毫米），菲涅尔近似的采样条件根本不成立。所以现在主流的全息再现实现里，算法库普遍用卷积法或角谱法代替直接菲涅尔变换。

3.2 卷积法实现：保持像素尺寸不变的核心技巧

卷积法守住了这么一条性质：把衍射积分写成冲激响应与物场的卷积，再利用卷积定理用两次 FFT 完成计算，输出平面像素尺寸和输入完全一致。全息再现的卷积形式对应传递函数

G(fx, fy) = exp[jkz·sqrt(1 - (λfx)² - (λfy)²)]

这个式子看着比菲涅尔积分更简洁，却保留了像面尺寸不变、任意 z 都能算的两个优势。代码里用scipy.signal.fftconvolve直接算，或者手写 FFT 相乘，两种方式等价。

from scipy.fft import fft2, ifft2, fftshift def angular_spectrum(u0, wavelength, pixel_size, z): """ 角谱法全息再现：像素尺寸不随 z 改变的衍射传播器。 参数 ---- u0 : 2D ndarray，传感器平面复振幅（频谱滤波后） wavelength : float，激光波长，单位米 pixel_size : float，传感器像素间距，单位米 z : float，光源/样品到传感器的再现距离，单位米 """ rows, cols = u0.shape fx = np.fft.fftfreq(cols, d=pixel_size) fy = np.fft.fftfreq(rows, d=pixel_size) FX, FY = np.meshgrid(fx, fy) # 频域传递函数：大于 1/λ 的部分 evanescent 波直接置 0 term = 1 - (wavelength * FX) ** 2 - (wavelength * FY) ** 2 term[term < 0] = 0 H = np.exp(1j * 2 * np.pi * z / wavelength * np.sqrt(term)) U = fft2(u0) * H result = ifft2(U) return result # 典型参数：氮氖激光 632.8nm，传感器像素 3.45μm，再现距离 85mm reconstructed = angular_spectrum(complex_field, wavelength=632.8e-9, pixel_size=3.45e-6, z=0.085)

参数选型上有三个坑值得展开。第一，wavelength * max(fx)若接近 1，则sqrt(term)的斜率很大，数值上对 z 极其敏感，微小的 z 误差会被放大成明显的离焦模糊；所以 z 的给定精度至少到毫米级，最好用自动对焦算法去搜。第二，网格点fx的排布是[-fs/2, fs/2)，零频在正中间，而ifft2默认零频在左上角，所以注意是否需要ifftshift，不匹配会导致输出像发生半周期平移。第三，当 z 很小而图像尺寸很大时，角谱法的频域采样间隔 δf=1/(N·Δ) 可能大于 1/λ 的两倍，此时存在混叠，请把原始全息图先裁到目标感兴趣区域再传播，而不是全幅传播。

3.3 再现像面参数的计算：放大率、像距和最小分辨尺寸

常见实验配置是物光经过显微物镜放大后再与参考光干涉，这时再现距离 z 与物镜焦距、筒长、管镜焦距有关。全息再现里最常用的标定公式是

M = z2 / z1

其中 z1 是物镜前焦面到样品的距离（近似等于物镜焦距），z2 是等效像距。数字全息系统的横向分辨率不取决于物镜数值孔径吗？既对也不对：NA 决定能收集多少衍射级次，但像素尺寸决定了频带能不能完整记录下来。Δx_sample = λ·z / (N·Δx_sensor)给出的是无显微放大下的样品面极限分辨尺寸。写代码时建议把上面这个公式保留在注释里，因为一旦遇到"为什么我重建的微球直径偏大 5%"这类问题，第一反应应该是去算系统放大率而不是去调滤波器半径。

4. 相位提取与解包裹：2π 跳变才是全息再现的最后一关

4.1 为什么距离信息藏在相位里，却读不出来

再现得到的复数场result每一项是 a + bj，相位直接np.angle(result)就行。但 arctan 的输出被截断在 (-π, π] 之间，而真实光程差往往大于一个波长，导致相位图上一圈一圈的条纹，相邻条纹之间差 2π。数字全息再现之后拿到的是包裹相位（wrapped phase），要变成连续的高度分布，必须做相位解包裹（phase unwrapping）。一维信号可以沿路径逐点累加 2π 偏移，二维图像却不行——噪声点会给路径积分引入 2π 整数倍的全局误差，这就是二维解包裹算法存在的理由。

4.2 质量图引导的路径解包裹：抗噪最稳的通用实现

最简单的逐行解包裹在条纹密集区（比如微球边缘）会整行漂移，产生胶皮状的伪影。我一般优先用质量图引导的洪水填充算法，原理是先计算每个像素的相位导数方差作为可靠性指标，从最可靠的像素出发，按队列依次展开邻域。

def phase_unwrap_quality(phase, reliability_mask=None): """ 基于质量图(二阶差分)的路径跟踪解包裹。 phase : 2D ndarray，包裹相位，范围(-π, π] reliability_mask : 可选，0~1 的权重图，用于屏蔽坏区域 """ # 相位导数的二阶差分 = 质量的负指标 dx = np.diff(phase, axis=1) dy = np.diff(phase, axis=0) # 处理 2π 跳变后求梯度 dx = np.angle(np.exp(1j * dx)) dy = np.angle(np.exp(1j * dy)) quality = -(np.abs(np.diff(dx, axis=0))[1:, :] + np.abs(np.diff(dy, axis=1))[:, 1:]) # 堆的起点是质量最高的像素；逐点弹出并扩展邻域 import heapq rows, cols = phase.shape unwrapped = np.zeros_like(phase) visited = np.zeros_like(phase, dtype=bool) heap = [] cy, cx = np.unravel_index(np.argmax(quality), quality.shape) # 第一行、最后一行、第一列、最后一列的像素质量不可靠，跳过做质量图的逻辑里处理 visited[cy, cx] = True unwrapped[cy, cx] = phase[cy, cx] heapq.heappush(heap, (0, cy, cx)) while heap: _, y, x = heapq.heappop(heap) for dy_, dx_ in ((1, 0), (-1, 0), (0, 1), (0, -1)): ny, nx = y + dy_, x + dx_ if 0 <= ny < rows and 0 <= nx < cols and not visited[ny, nx]: delta = phase[ny, nx] - phase[y, x] delta -= 2 * np.pi * np.round(delta / (2 * np.pi)) unwrapped[ny, nx] = unwrapped[y, x] + delta visited[ny, nx] = True # 质量图的索引比相位图小1，要小心对齐 q = quality[ny - 1, nx - 1] if ny > 0 and nx > 0 else -999 heapq.heappush(heap, (-q, ny, nx)) return unwrapped unwrap_phase = phase_unwrap_quality(np.angle(reconstructed))

这段代码里的delta -= 2π * round(delta / 2π)是核心：它保证两个像素之间的相位差被调整到最短路径上，无论相位图上这一格是 3.1 还是 -3.1，相对关系都不会偏。质量图用二阶差分反映局部条纹密度，密度越大质量越低，扩张就绕过这些区域。参数quality与phase的尺寸差一像素，代码里用if ny > 0 and nx > 0兜底，但如果输入图像较大，更可靠的做法是把质量图先np.pad到与相位图同尺寸。

4.3 从解包裹相位到高度：载波相位的扣除方法和标定公式

全息再现得到相位差之后，还不能直接乘系数。干涉图里的条纹还包含参考光和物光夹角带来的线性载波相位（tilt phase），这一项不扣掉，解出来的高度会叠加一个斜坡面。常见的扣除办法是拿无样品的纯平场全息图做同样的再现与解包裹，得到载波相位分布，再用样品的解包裹相位减去它。

# 假设 flat_unwrapped 是平场全息图用同样流程解包裹后的相位 tilt = flat_unwrapped # 线性斜坡，但含低频偏移 corrected_phase = unwrap_phase - tilt # 高度换算：氮氖激光反射式测量，光程差为 2h wavelength = 632.8e-9 height = corrected_phase * wavelength / (4 * np.pi)

这里有坑：透射式还是反射式，换算系数差一倍。透射式测量样品折射率未知时分不清是厚度还是折射率变化，反射式最干净，公式里除以 4π 因为光走了来回。corrected_phase要做趋势去除再去算高度，但趋势去除时不要用多项式拟合整个面——样品占全场 30% 以下时，应该只选四角背景区域来拟合平面，否则把样品的真实弯曲一起"拟合"掉了。

提示：解包裹相位里出现孤立的±2π斑块，通常不是物理信号，而是传感器坏点或掩模边缘的频谱截断伪影。在质量图里把这些像素权重设 0，比解包裹后再做中值滤波更有效。

5. 数字全息再现的验证手段与现场检查技巧

5.1 三步验证：平面镜标定、已知高度台阶、重复性

拿到一个全息再现系统，第一步不是测量未知样品，而是用光学平晶或镀膜硅片当样品记录。理想平面在解包裹后的相位图上应该是平坦的，残余起伏就是系统像差。记录平面镜全息图，做完整链路（滤波、角谱传播、解包裹），输出残余相位峰谷值；数字全息系统的可用横向分辨率就是峰谷值对应高度的两到三倍。第二步用刻蚀台阶或标准高度块（高度已知到纳米级）验证高度比例系数；这时就会暴露波长校准问题，激光器标称 632.8nm 但实际可能是 632.99nm，校准系数直接乘到最终高度上。

第三步是重复性，同一位置连续拍十张，解包裹后计算每张相位的标准差。这个值告诉你系统的测量噪声是 0.5nm 还是 5nm，也直接决定你论文里误差棒该画多宽。这里有一项值得做：把十张全息图先平均再处理，和平场平均一样，系统随机噪声会按 1/√N 下降。

5.2 离焦判据：用锐度曲线而不是肉眼

自动对焦搜索 z 时，很多新手用强度图最锐作为判据，但相位图对焦更敏感。我常用的做法是定义相位梯度能量，或直接用高频傅里叶能量作为对焦评价函数。具体命令如下

def focus_metric(recon): grad = np.gradient(np.angle(recon)) return np.sum(grad[0] ** 2 + grad[1] ** 2) z_best = None score_best = -1 for z in np.linspace(0.080, 0.090, 50): r = angular_spectrum(complex_field, 632.8e-9, 3.45e-6, z) s = focus_metric(r) if s > score_best: score_best, z_best = s, z

np.gradient在 10 像素边缘会放大噪声，所以曲线在近焦点附近会变尖；扫描步长取 0.2mm 足够。值得补充一点：用相位梯度作测度时，样品边缘的相位跳变贡献很多评分，如果样品本身是阶梯状，焦点位置会自动偏向使阶梯最陡的地方——这时要和幅度图的对比度做交叉验证。

5.3 快速故障排查表：当你"什么都重建不出来"

最后一条经验，把最常见的失败现象与对应处置集中在一起，方便现场翻。

磁盘上写完这么多，记得把每一步处理后的复数场和数据参数同时保存成.npz，而不是只存最终 PNG。因为调整滤波半径或 z 参数的代价是重新跑一遍后端，有了中间结果，调试周期从两小时缩短到两分钟。这个习惯比任何算法技巧都更能决定全息再现项目能不能按期交付。

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