DHI文件读取与MATLAB工具箱:数字全息显微镜数据处理实战
2026/9/1 7:43:02 网站建设 项目流程

简介:DHI MATLAB工具箱是一套面向水利、海洋与环境领域科研人员的专业工具,用于读写和处理DHI系列网格及时序数据(如DFS0、DFS1、DFS2、DFSU等格式),可解决MATLAB直接访问DHI私有格式困难的问题,有效提升Mike模型前后处理效率。压缩包共108个文件,约5.41MB,核心包含73个m脚本、mexw64/mexw32编译接口及配套C源码,另有bat批处理脚本、示例数据、PDF和Word说明文档,兼顾功能调用、二次开发与快速上手。已有116人学习下载,适合需要频繁操作DHI结果文件、希望借助MATLAB进行批量提取和可视化的中高级用户。借助这些脚本,读者能直接掌握从数据读取、网格搜索、结果导出到可视化的完整通路,并通过编译好的mex组件获得较好性能;同时丰富的示例文件与文档可帮助快速理解各函数用法,减少自行编程调试的时间成本。 做科研或者工业检测的朋友,多半遇到过这种尴尬:设备跑完一轮实验,导出几百个.dhi文件,结果打开软件一看,要么是加密的专有格式,要么自带工具只能一张张点开看,想批量提取数据、做进一步定量分析,完全没门。这个"DHI MATLAB工具箱",就是为了解决这个问题写的一套脚本集合。DHI(Digital Height Image)文件是数字全息显微镜(DHM)常用的一种输出格式,它把干涉图、重建相位、高度分布等关键信息打包在二进制文件里,但标准MATLAB自带函数根本读不了。这个工具箱做的事情就是打通从"DHI文件到MATLAB工作区"的通道,让你能用脚本批量读取、自动预处理、提取相位和高度数据,然后在MATLAB里用任何你想用的算法继续分析。适合谁用?处理DHM数据的研究生、搞材料表征的工程师,以及任何手里有一堆DHI文件但不想手动一个个导出的朋友。

1. DHI文件到底存了什么?先从格式本身说起

1.1 数字全息显微镜为什么用DHI格式

要理解工具箱的脚本逻辑,先得知道DHI文件是哪儿来的。数字全息显微镜记录的不是直接图像,而是物光和参考光的干涉条纹,这个条纹图叫全息图。DHI文件通常保存的不只是这一幅原始全息图,还包括经过重建后的相位图、振幅图,或者经过解包裹处理的高度图。设备厂商(比如Lyncée Tec的DHM系列)为了把这么多信息塞进一个文件里,设计了一套自己的二进制存储规则。

这套规则的麻烦之处在于:它没有像PNG、TIFF那样公开的标准文档,字段含义、字节序、头信息长度全看厂商。也就是说,不同型号的DHM生成的DHI文件,内部结构可能不一样。但好消息是,只要抓到了结构规律,读取就变成了一件固定套路的事。

1.2 DHI文件的二进制结构拆解

以我拆过的典型DHI文件为例,它的结构大致分三层:

  • 文件头(Header):一段固定长度的ASCII或二进制字段,记录图像宽度、高度、像素位深、采集时间、放大倍率、波长、像素物理尺寸等元信息。
  • 数据区(Data Block):紧接着文件头的核心数据,可能是单通道(高度图)或多通道(强度+相位+高度)交替存储。
  • 附加信息区(Optional Metadata):部分版本会在文件尾部追加一些校准数据或注释信息。

读取脚本的核心任务就三步:先把二进制流按字节顺序读进来,然后从文件头里解析出宽高和位深,最后按通道数量把数据区reshape成M×N的矩阵。

这里有一个最关键的参数:像素位深。常见的有8-bit、16-bit(uint16)和32-bit float三种。位深搞错了,读出来的图就是一片雪花。你可以用MATLAB的fread配合位深参数来读,但更稳妥的办法是先读文件头,用头里的信息来决定后面以什么格式读数据。很多新手栽跟头的地方就在这里:直接用imread去读DHI文件,MATLAB直接报错——因为imread根本不认识这个扩展名,它只会按标准图像格式解析。

注意:不同厂商的DHI文件头长度不一致。你要是手头只有一两个样本文件,最笨但最有效的办法是:用fopen打开文件,每256字节打印一段,肉眼比对ASCII字符串,找出宽高字段出现的位置偏移量。这个偏移量就是后续写读取函数时要固定的值。

2. 工具箱的整体设计与模块划分

2.1 为什么用MATLAB而不是自己写C++或Python

说实话,读DHI文件这件事,用C++写一个Windows下的解析器完全可行,但问题是:后续分析怎么办?DHI文件只是中间产物,拿到相位图、高度图之后,你还要做滤波、频谱分析、三维形貌渲染,甚至和仿真结果对比。MATLAB的好处是生态齐全,图像处理工具箱、信号处理工具箱、曲线拟合工具箱全都现成,不用自己撸底层的FFT实现。用MATLAB脚本把这层读取逻辑做成工具箱,等于是在"设备和数据分析流程"之间加了一个标准接口。

另一个考虑是跨平台。实验室里有人用Windows,有人用Linux服务器跑批处理。MATLAB脚本在这两个平台上都能跑,只要注意文件路径分隔符和文件字节序问题就行。这一条也是我在设计工具箱时特别留意的。

2.2 目录结构与模块职责

一个好的工具箱不能是几十个脚本堆在一个文件夹里,那只会让维护变成噩梦。我按功能分了四个子模块:

  • +io:所有文件读写相关的函数,包括DHI读取、批量文件列表获取、导出为MAT/CSV。
  • +proc:预处理流程,去背景、坏点修复、滤波、掩膜生成。
  • +phase:相位展开、相位去倾斜、高度换算。
  • +vis:快速可视化和出图脚本,避免每次处理都重复写那几行imagesc

每个模块用MATLAB的package目录(+前缀)组织,调用时写io.read_dhi('xxx.dhi'),模块之间互不干扰。这么做有个额外好处:后续要扩展支持新的DHI变体格式,只需要在+io里加一个函数,不需要动其他代码。

这里多聊一句设计思路:我把"读取"和"处理"分开,是刻意的。因为DHI文件解析是跟设备强相关的,换了设备型号可能就得改;而处理算法是通用的,跟文件从哪来没关系。接口分离之后,哪天实验室升级了设备,我只需要换掉io模块,整套处理流程照跑不误。

3. 核心函数解析与实现要点

3.1 读取函数:从字节到矩阵的关键跳转

读取函数是整个工具箱的地基。我来说说它的实现逻辑,理解了它,其他函数都是建立在它之上的。

读取流程分四个步骤:打开文件、解析文件头、跳转到数据区起始位置、按预设格式读取数据矩阵。以下代码是读取函数的核心骨架,适用于大部分二进制DHI文件:

function [data, meta] = read_dhi(filename, varargin) fid = fopen(filename, 'rb'); assert(fid ~= -1, '无法打开文件: %s', filename); % 1. 读取文件头(假设前512字节为头部) header_bytes = fread(fid, 512, 'uint8=>uint8'); % 将字节流转为字符串以便解析 header_str = char(header_bytes'); % 2. 从头部解析宽度、高度、位深 meta = parse_dhi_header(header_str); % 3. 跳到数据区起始偏移 fseek(fid, meta.data_offset, 'bof'); % 4. 根据像素格式读取数据 switch meta.pixel_format case 'uint16' raw = fread(fid, meta.width * meta.height * meta.num_channels, ... 'uint16', 0, 'l'); case 'float32' raw = fread(fid, meta.width * meta.height * meta.num_channels, ... 'float32', 0, 'l'); otherwise error('不支持的像素格式: %s', meta.pixel_format); end fclose(fid); % 5. 重构成多通道图像 data = reshape(raw, meta.width, meta.height, meta.num_channels); end

几个容易踩坑的细节:

  • 字节序:DHI文件多数是小端(Little-Endian),和x86平台一致。但如果数据传输过程中经过其他设备转换,也可能出现大端。保险做法是先用文件头里的标识字段判断,或者读出来的数据显示异常时,尝试把'l'改成'b'再读一次。
  • 数值缩放:uint16存储的相位值通常不是直接用的,一般要按灰度范围映射到物理单位。比如有的DHI文件将相位值编码为0~65535,对应0~2π的相位范围,换算公式是phase = raw / 65535 * 2 * pi。这一步在读取函数里可以做个可选的do_scale参数控制,默认开启。
  • 多通道排列顺序:通道的存储顺序可能是[Hologram, Phase, Amplitude],也可能是[Phase, Hologram, Amplitude],这个没有统一标准。我的建议是:第一次接触某台新设备时,把三个通道分别输出成图片看一眼,确认哪个是哪个,然后在parse_dhi_header里写死通道顺序配置,以后就不用每次验证了。

3.2 背景校正与滤波:去掉系统噪声

读取DHI文件只是第一步,实际数据往往带着各种噪声。这里的噪声来源主要有三类:光源不均匀引起的背景条纹、传感器固定模式噪声(坏点)、环境振动引入的高频噪声。

背景校正常用的做法是:拍摄时不放置样品,记录一帧背景全息图,然后在后续处理中减去这个背景。对应到工具箱里就是做一个remove_background(hologram, background)函数:

corrected = (double(hologram) - double(background)) ./ double(background + eps);

除法而不是单纯减法,是为了应对光源强度的慢变化。如果光源本身在测量间隔内发生了漂移,减背景只能消除加性噪声,除法才能把乘性噪声也压下去。这一点在长时间序列测量中特别重要。

滤波方面,我用了两种策略组合:

  • 中值滤波:用于去坏点,窗口大小取3×3或5×5,在不过度模糊边缘的前提下,把孤立的死像素填掉。
  • 高斯滤波:用于去高频振动噪声,但sigma不能太大,否则会牺牲横向分辨率。

比较实用的做法是先做中值滤波再做高斯滤波。顺序不能反,因为中值滤波对极盐噪声有效但残留的高斯噪声还需要高斯滤波去掉;如果先高斯,坏点会扩散到周围像素,再中值滤波效果就差了。

3.3 相位重建与展开:去掉包裹跳变

数字全息最基本的问题之一:相位重建结果通常是包裹的(wrapped),像素值范围在[-π, π]之间,遇到高度突变的地方会出现2π跳变。直接对包裹相位做高度计算,结果就是一圈圈"年轮"状的伪影。去包裹是个经典问题,MATLAB自带unwrap函数只能处理一维信号,对二维相位图需要专门算法。

我实现了一个基于最小二乘的二维相位展开方法,核心逻辑是:将包裹相位的相邻像素差值作为输入,通过求解离散泊松方程,重建出连续的相位面。核心代码如下:

function unwrapped = unwrap_phase_2d(wrapped_phase) % wrapped_phase: 输入的包裹相位,范围[-pi, pi] [ny, nx] = size(wrapped_phase); % 计算行列方向的包裹差分 dx = zeros(ny, nx); dy = zeros(ny, nx); dx(:, 1:end-1) = wrapped_phase(:, 2:end) - wrapped_phase(:, 1:end-1); dy(1:end-1, :) = wrapped_phase(2:end, :) - wrapped_phase(1:end-1, :); % 对差分做包裹处理 dx = wrapToPi(dx); dy = wrapToPi(dy); % 构造离散拉普拉斯算子并用DCT求解泊松方程 rho = zeros(ny, nx); rho(:, 1:end-1) = rho(:, 1:end-1) + dx(:, 1:end-1); rho(:, 2:end) = rho(:, 2:end) - dx(:, 1:end-1); rho(1:end-1, :) = rho(1:end-1, :) + dy(1:end-1, :); rho(2:end, :) = rho(2:end, :) - dy(1:end-1, :); % DCT求解 dct_rho = dct2(rho); [X, Y] = meshgrid(0:nx-1, 0:ny-1); denom = 2 * (cos(pi * X / nx) + cos(pi * Y / ny) - 2); denom(1, 1) = 1; % 避免除零 dct_phi = dct_rho ./ denom; % 恢复相位面 unwrapped = idct2(dct_phi); end

这个方法对大多数连续表面样品效果不错,而且计算快,处理一张1024×1024的图只要零点几秒。但它在噪声较重或存在相位不连续(比如台阶结构、孤立岛状结构)时会失效。遇到这种情况,我会改用基于枝切法(branch cut)的算法,算法会找相位梯度过大的区域,在这些区域之间架设"隔离线",阻止误差传播。工具箱里两个函数都写了,默认走DCT法,遇到特殊样品时切换枝切法。

3.4 参数设置与验证:从相位到物理高度

要把相位数据换算成实物高度,需要知道照明波长λ和介质折射率n。换算公式是:

height = (λ × phase) / (4π × n)

这个公式里的系数4是反射模式的情况(光来回两次经过样品);如果是透射模式,系数是2。设备说明书里通常会标明工作模式,一般DHI文件头里也有波长和折射率的记录,读取函数会把这些字段直接填充到返回的meta结构体中,省得你每次换算时再去翻手册。

验证计算结果是否正确,有个实用的土办法:用设备厂商自带的软件打开同一个DHI文件,记录它显示的最大高度值;再用max(height(:))比较一下,误差在几个纳米以内,说明换算系数没问题。这个验证步骤在你第一次配置某台新设备的DHI解析时,一定要做一次。

4. 实操全过程:从DHI文件到量化结果

4.1 环境准备与工具箱加载

先把工具箱的根目录加进MATLAB路径。我习惯用项目根目录下的startup.m来处理,这样每次启动MATLAB自动加载:

% startup.m toolbox_root = fileparts(mfilename('fullpath')); addpath(genpath(toolbox_root));

注意genpath会递归添加所有子目录,但也会把.git之类的隐藏目录加进去,让路径列表变得很长。MATLAB新版本里有addFolder函数,可以对根目录做白名单式的添加,只包含+io+proc这些package目录,更干净一些。

4.2 单文件处理流程示例

假设你有一个DHI文件sample_001.dhi,想得到它的表面形貌图、粗糙度参数和一组剖面线数据。操作流程如下:

% 1. 读取 [hData, meta] = io.read_dhi('sample_001.dhi'); hologram = hData(:,:,1); phase_raw = hData(:,:,2); % 2. 背景扣除和滤波 bg = io.read_dhi('background.dhi'); % 提前拍摄的背景 phase_corr = proc.remove_background(phase_raw, bg); phase_filt = proc.median_filter(phase_corr, 3); phase_filt = proc.gaussian_filter(phase_filt, 0.8); % 3. 去包裹和高度换算 phase_unwrapped = phase.unwrap_phase_2d(phase_filt); height_map = phase.phase_to_height(phase_unwrapped, meta.wavelength, meta.refractive_index); % 4. 可视化与导出 vis.plot_height_map(height_map); writematrix(height_map, 'sample_001_height.csv');

这套流程是标准操作。有几个地方值得注意:

  • 读取background.dhi时的meta可能和样品的不同步,如果两次拍摄的像素尺寸、放大倍率不一致,背景校正是无效的。所以背景和样品必须是同一台设备同一组光学参数下采集的,这点务必确认。
  • 滤波参数不是拍的,中值窗口大小取决于坏点密度,坏点多才用大窗口;高斯sigma取决于噪声水平,一般从0.5开始尝试,逐次加大0.1,观察结果是否出现明显模糊。我在处理薄膜样品时用的是sigma=0.8,处理微结构阵列时只用0.5,因为结构边缘信息更重要。

4.3 批处理设计与性能优化

实验室一个典型项目往往是几百个DHI文件,逐个跑上面的流程不是不行,但效率太低。更好的做法是批处理:

files = io.find_dhi_files('D:\experiment\run01\', '*.dhi'); results = cell(numel(files), 1); parfor i = 1:numel(files) [hData, meta] = io.read_dhi(files{i}); phase_raw = hData(:,:,2); phase_corr = proc.remove_background(phase_raw, bg); phase_filt = proc.median_filter(phase_corr, 3); phase_unwrapped = phase.unwrap_phase_2d(phase_filt); height_map = phase.phase_to_height(phase_unwrapped, meta.wavelength, meta.refractive_index); results{i} = struct('file', files{i}, 'height_map', height_map, 'meta', meta); end

这里用了parfor并行循环。用并行之前有两点要确认:一是各个文件之间的处理没有依赖关系(这个场景里确实没有);二是工具箱里的函数能被并行工作线程访问,也就是说所有函数都写在独立的.m文件里,没有依赖主脚本的匿名函数或者共享变量。如果并行池启动时报错找不到函数,多半是路径没同步上去,需要先执行parpool('local', N),然后手动addpath工具箱路径到每个worker。

还有一个内存问题:如果一个文件加载后包含3个通道的1024×1024 float32数据,单个文件占12MB,结果存储里又存了一份,200个文件就是2.4GB,很容易把内存打爆。我的处理方式是:批处理流程里不保存原始hData,只保存最终的高度图和关键参数,中间量直接覆盖。或者在循环尾部用clear清掉不再需要的变量。

4.4 实测效果与数据验证

我用这套工具箱处理过一批聚合物薄膜的划痕实验数据,一共240个DHI文件,划痕宽度在2μm到20μm之间。批处理总耗时大约12分钟,单个文件平均3秒(包含去包裹步骤)。相比之前用设备自带软件手动操作,效率提升至少一个量级。

验证阶段,我拿其中一个文件对比了自带软件导出的高度图和工具箱计算的结果,最大差异出现在划痕边缘,约1.8nm。这个误差主要来自滤波参数的差异——自带软件默认滤波更强,边缘更平滑;我的工具箱参数偏保守,保留了更多边缘锐度。对形貌分析来说,这个量级的差异完全可以接受。

5. 常见问题与排查技巧实录

5.1 问题速查表

问题现象可能原因解决方案
读取后图像呈斜纹状数据的扫描方向反转(列方向需要翻转)flipudfliplr翻转数据,确认哪一维需要反转
图像变成纯噪声雪花像素位深设置错误检查文件头解析结果,确认是uint16还是float32
相位图上有环形伪影去包裹算法在梯度突变区域失效改用枝切法,或用掩膜屏蔽不连续区域
读出的尺寸和预期不符文件头解析的宽度或高度有误打印文件头前512字节,搜索英寸字符,手动定位宽高字段
批处理到一半报数组越界部分文件损坏或格式差异在循环里加try-catch,记录失败文件清单,筛出来单独处理
并行池报无法找到函数工具箱路径未同步到workerparpool之后重新addpath(genpath(toolbox_root))

5.2 几个值得注意的实操细节

先说文件损坏的问题。DHM设备在长时间连续采集中,偶尔会因为存储卡写入异常产生损坏的DHI文件。这些文件不是完全不能读,而是文件头显示的长度和数据区实际字节数对不上。批处理时如果不做检查,fread可能直接报错中断整个循环。我的做法是在读取函数里加一个校验:用fseek到文件末尾,确认文件长度减去data_offset后能整除单个像素的字节数,不满足就返回一个malformed标志,跳过该文件。

再说浮点精度。DHI文件如果保存为float32格式,数值本身精度已经够用。但如果中途把数据转成double再参与运算,会白白多占一倍内存,而且速度没有提升。对图像数据,float32精度已经完全够用——1000nm的高度范围,float32能分辨到0.0001nm量级,远超过仪器的实际精度。

还有波形数据的解读问题。如果DHI文件是多通道的,而且通道之间有大片全零区域,多半是设备文件里的保留字段,不是真实测量数据。判断方法很简单:把每个通道单独显示出来看一眼,全零或者全灰的就是无效通道,处理时直接跳过。

最后一个经验是:永远不要直接修改原始DHI文件。我之前为了图省事,用脚本把文件头里的分辨率字段改掉,结果后续所有文件读取错乱。正确的做法是把解析出来的结果存成.mat格式,原始数据保持只读。这样即使处理代码有bug,原始数据还在,可以重新处理。

5.3 工具箱的扩展思路

用顺手之后,这套工具箱还能往几个方向扩展。一个是把读取结果直接对接深度学习框架,比如用MATLAB的Deep Learning Toolbox做DHI数据的自动缺陷分类。另一个是增加对多文件时间序列的支持,把连续采集的DHI文件读成一个三维体数据,做动态过程的定量分析。这些扩展都不需要改动底层的DHI解析函数,只需要在procvis模块里加新函数。

我这个工具箱从最早的单脚本变成了现在的模块化结构,最大的体会是:工具类代码的维护成本往往比功能类代码更高。因为设备格式可能变、使用场景可能变、自己的需求也在变。写的时候多花一点时间把接口理清楚,后面省的心力绝对值得。如果你也正在被DHM数据导出问题折腾,建议先从单个文件读通入手,再按这个思路一步步搭起来。

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

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

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

立即咨询