简介这份资源是一套基于MATLAB实现的数字图像相关DIC主程序面向材料科学、工程力学与生物医学领域的学习者和研究人员用于从连续图像序列中非接触式地计算物体表面的位移、应变与应力。压缩包共161个文件包含82个m脚本、43张jpg与19张png图像、10个html说明页、5个fig图形文件及少量txt文档整体约6.01MB代码、示例图与界面文件配套齐全。资源围绕图像预处理、特征点检测与配对、图像匹配、位移场计算、应变应力求解以及结果可视化等环节展开脚本与函数分工明确便于读者按模块理解DIC算法的完整流程。目前已有4164人学习下载适合希望借助MATLAB上手DIC、调试算法或将其迁移到自身实验场景的读者参考与二次开发。1. 从一组拉伸试验照片说起这套 DIC 主程序到底能算什么手里有一组材料拉伸试验的连拍照片想拿到全场位移、应变甚至应力分布传统做法是贴应变片、逐点测量费时费力还只能拿到离散几个点的数据。这套 DICDigital Image Correlation数字图像相关主程序就是干这个的输入一叠按时间顺序排列的图像序列输出每一帧相对于参考帧的位移场、由位移场差分得到的应变场再结合材料本构关系换算出应力分布。它面向的是做材料力学实验、结构变形监测、疲劳裂纹观测这类场景的开发者尤其是习惯用 MATLAB 做数据处理和二次开发的人。整套代码是纯 MATLAB 实现不依赖额外工具箱也能跑基础流程适合拿来当 DIC 算法学习模板也适合直接改成自己实验台的批处理工具。2. 环境准备与目录结构把主程序跑起来前要理清的东西2.1 运行环境与依赖清单这套程序对 MATLAB 版本要求不算苛刻R2018b 及以上基本都能跑核心用到的函数集中在图像处理工具箱Image Processing Toolbox里比如imread、imresize、normxcorr2、imgaussfilt这些。如果只做最基础的整像素相关理论上不装工具箱也能用conv2手写互相关但速度会慢一个量级常见做法还是老老实实把工具箱装上。内存方面要注意DIC 的搜索窗口和子区尺寸直接决定内存占用。假设图像是 2048×2048子区 41×41搜索半径 20 像素单帧相关计算会产生一个 2048×2048×41×41 量级的中间矩阵如果一次性把所有帧读进内存8GB 内存的机器很容易被撑爆。我一般会强制逐帧读取、逐帧释放用clear和pack控制内存峰值。% 环境自检脚本 check_env.m assert(~isempty(ver(images)), 缺少图像处理工具箱请先安装); img imread(frame_0001.tif); assert(size(img,3)1, 当前版本仅支持灰度图彩色图请先转灰度); fprintf(图像尺寸: %dx%d, 数据类型: %s\n, size(img,1), size(img,2), class(img));这段自检脚本做了三件事确认工具箱存在、确认输入是单通道灰度图、打印图像基本信息。参数上frame_0001.tif是参考帧文件名实际使用时替换成自己序列的第一帧即可。如果读进来是彩色图后面相关计算会按三通道分别处理结果完全错乱所以这一步断言不能省。2.2 目录结构与文件职责拿到源码包后先别急着点运行把目录结构过一遍能省掉后面大量找文件的时间。典型的组织方式是这样DIC_Main/ ├── main.m % 主入口串起整个流程 ├── config/ │ └── params.m % 所有可调参数集中在这里 ├── core/ │ ├── dic_corr.m % 互相关核心算整像素位移 │ ├── subpixel.m % 亚像素拟合精度到 0.01 像素 │ ├── strain_calc.m % 位移场差分求应变 │ └── stress_calc.m % 本构关系换算应力 ├── io/ │ ├── load_sequence.m % 批量读图 │ └── export_result.m % 结果导出为 mat/csv ├── utils/ │ ├── roi_mask.m % 感兴趣区域掩膜 │ └── plot_field.m % 位移/应变云图绘制 └── data/ └── sample/ % 示例图像序列main.m是唯一需要手动改路径的地方其余模块通过函数调用串联。params.m把所有魔法数字集中管理改子区大小、搜索半径、应变窗口这些参数不用翻遍每个文件。这种分层的好处是调算法只动core/换数据格式只动io/画图不满意只动utils/互不干扰。2.3 参数配置几个决定成败的数字params.m里有几个参数直接决定结果好坏新手最容易在这里翻车参数名含义典型值影响subset_size子区边长像素31~51太小噪声大太大空间分辨率低search_radius搜索半径像素15~30必须大于预估最大位移step_size计算点间隔5~10越小点越密计算量平方增长strain_window应变差分窗口5~11太小应变噪声大太大平滑过度pixel_size物理像素尺寸mm由标定得决定位移物理量换算% params.m 片段 params.subset_size 41; % 子区 41x41兼顾噪声和分辨率 params.search_radius 25; % 搜索半径 25 像素覆盖预估位移 params.step_size 5; % 每 5 像素算一个点 params.strain_window 7; % 7 点差分求应变 params.pixel_size 0.02; % 每像素 0.02 mm需标定得到 params.roi [200 200 1800 1800]; % [x y w h] 感兴趣区域search_radius是最容易设错的一个。如果材料实际位移超过搜索半径相关峰会跑到搜索窗口外面程序要么报错要么给出完全错误的位移。稳妥做法是先跑一帧看最大位移大概多少再把搜索半径设成它的 1.5 倍。pixel_size必须通过标定板或已知尺寸参照物得到直接填 1 会让位移单位变成像素而不是毫米后面应力计算全错。3. 核心算法拆解从互相关到应力场的完整链路3.1 整像素互相关dic_corr.m在算什么DIC 的起点是模板匹配。取参考帧上一个子区比如 41×41 的方块在变形帧的搜索窗口里滑动找归一化互相关NCC最大的位置那个位置就是子区中心的整像素位移。核心代码大致长这样function [u, v, cc_max] dic_corr(ref_img, def_img, params) % 输入: 参考帧、变形帧、参数结构体 % 输出: 整像素位移 u(水平) v(垂直)、相关系数峰值 half floor(params.subset_size / 2); [H, W] size(ref_img); % 计算点网格 xs params.roi(1)half : params.step_size : params.roi(1)params.roi(3)-half; ys params.roi(2)half : params.step_size : params.roi(2)params.roi(4)-half; u zeros(numel(ys), numel(xs)); v zeros(numel(ys), numel(xs)); cc_max zeros(numel(ys), numel(xs)); for iy 1:numel(ys) for ix 1:numel(xs) % 提取参考子区 r ref_img(ys(iy)-half:ys(iy)half, xs(ix)-half:xs(ix)half); % 在搜索窗口内做归一化互相关 sr params.search_radius; win def_img(ys(iy)-half-sr:ys(iy)halfsr, ... xs(ix)-half-sr:xs(ix)halfsr); cc normxcorr2(r, win); [cc_max(iy,ix), idx] max(cc(:)); [py, px] ind2sub(size(cc), idx); % 换算成位移注意 normxcorr2 输出尺寸偏移 u(iy,ix) px - (size(cc,2) - half - sr) - 1; v(iy,ix) py - (size(cc,1) - half - sr) - 1; end end end逐段说明half是子区半宽xs/ys是按step_size生成的网格点坐标只保留子区完全落在 ROI 内的点。双重循环里r是参考子区win是变形帧上以当前点为中心、向外扩search_radius的搜索窗口。normxcorr2返回的相关系数矩阵尺寸是(subset2*sr) × (subset2*sr)峰值位置(py,px)需要减去窗口中心偏移才能得到真实位移。这个偏移量算错是新手最常见的 bug表现为位移整体偏了一个常数。cc_max是质量指标一般要求大于 0.8低于 0.6 的点基本不可信后面应变计算时应该剔除或标记。实际项目中我会在export_result.m里把cc_max一起导出方便事后筛点。3.2 亚像素拟合把精度从 1 像素推到 0.01 像素整像素相关只能给到 1 像素精度对大多数材料试验来说远远不够。subpixel.m在整像素峰值附近做曲面拟合常用的是二次多项式拟合或高斯拟合。二次拟合实现简单、速度快是工程上最常见的选择function [du, dv] subpixel(cc, py, px) % 输入: 相关系数矩阵、整像素峰值行列 % 输出: 亚像素修正量 du dv % 边界保护 if py1 || pysize(cc,1) || px1 || pxsize(cc,2) du 0; dv 0; return; end % 水平方向三点二次拟合 c1 cc(py, px-1); c2 cc(py, px); c3 cc(py, px1); denom 2*(c1 - 2*c2 c3); if abs(denom) 1e-12 du 0; else du (c1 - c3) / denom; end % 垂直方向同理 r1 cc(py-1, px); r2 cc(py, px); r3 cc(py1, px); denom 2*(r1 - 2*r2 r3); if abs(denom) 1e-12 dv 0; else dv (r1 - r3) / denom; end % 限制修正量在 ±0.5 像素内防止拟合异常 du max(min(du, 0.5), -0.5); dv max(min(dv, 0.5), -0.5); end二次拟合的公式来自对三个相邻相关系数点求抛物线顶点du就是顶点相对中心点的偏移。denom接近零说明三点几乎共线拟合无意义直接返回零修正。最后把修正量钳制在 ±0.5 像素内是因为理论上亚像素修正不应该超过半个像素超出说明相关峰形状异常多半是噪声或遮挡导致的。这一步做完位移精度通常能到 0.01~0.05 像素配合合适的pixel_size物理位移精度可以到亚微米级。3.3 应变场计算差分窗口怎么选拿到每个网格点的位移后应变就是位移的空间导数。strain_calc.m用局部最小二乘或简单中心差分实现function [exx, eyy, exy] strain_calc(u, v, step, win) % 输入: 位移场 u v、网格间距 step(像素)、差分窗口 win(点数) % 输出: 正应变 exx eyy、剪应变 exy [ny, nx] size(u); exx zeros(ny, nx); eyy zeros(ny, nx); exy zeros(ny, nx); h win * step; % 物理差分距离 for i 1:ny for j 1:nx i1 max(1, i-win); i2 min(ny, iwin); j1 max(1, j-win); j2 min(nx, jwin); % 局部平面拟合求偏导 [dudx, dudy] local_grad(u, i1,i2,j1,j2, step); [dvdx, dvdy] local_grad(v, i1,i2,j1,j2, step); exx(i,j) dudx; eyy(i,j) dvdy; exy(i,j) 0.5*(dudy dvdx); end end endwin控制差分窗口大小窗口越大应变越平滑但空间分辨率越低。local_grad内部用最小二乘拟合一个平面比简单中心差分抗噪能力强。这里有个血泪经验应变场对位移噪声极其敏感位移误差 0.01 像素、网格间距 5 像素时应变噪声就有 0.002已经和很多材料的弹性应变一个量级。所以strain_window不能设太小通常 5~11 之间具体看噪声水平和关心的应变尺度。3.4 应力换算本构关系怎么接进来应力不是图像直接测出来的而是由应变通过材料本构关系算出来的。stress_calc.m支持线弹性和简单塑性两种模式function [sxx, syy, sxy] stress_calc(exx, eyy, exy, mat) % 输入: 应变场、材料参数结构体 % 输出: 应力场 if strcmp(mat.model, elastic) % 平面应力线弹性 E mat.E; nu mat.nu; C E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; sxx C(1,1)*exx C(1,2)*eyy; syy C(2,1)*exx C(2,2)*eyy; sxy C(3,3)*exy; elseif strcmp(mat.model, plastic) % 增量塑性需要加载历史 [sxx, syy, sxy] plastic_update(exx, eyy, exy, mat); end end线弹性部分就是标准的平面应力刚度矩阵乘应变向量E是弹性模量nu是泊松比。塑性模式需要传入加载历史因为塑性应变和加载路径有关不能只看当前应变。实际项目中如果只是做弹性阶段验证用线弹性模式就够了要跟踪屈服和硬化就得把整个加载序列的应变按时间顺序喂进去逐步更新。4. 避坑与常见问题跑不通时先查这几条4.1 位移整体偏一个常数现象所有点的位移都比预期大或小一个固定值。原因normxcorr2输出矩阵的坐标原点和子区中心的换算关系搞错了这是最经典的翻车点。解决用一张已知平移量的合成图验证比如把参考图整体右移 3 像素生成变形图看程序输出是不是 3。如果不是检查dic_corr.m里px - (size(cc,2) - half - sr) - 1这个偏移公式不同 MATLAB 版本normxcorr2的边界处理略有差异建议用合成图实测标定一次。4.2 相关系数普遍偏低现象cc_max大部分点低于 0.6位移场看起来像噪声。原因通常有三个一是子区太小纹理信息不足二是散斑质量差对比度低或颗粒太细三是搜索半径不够真实位移超出搜索范围。解决先把subset_size从 31 加到 51 试试如果还不行检查散斑图是不是过曝或欠曝必要时重新制斑。搜索半径按预估最大位移的 1.5 倍设置。4.3 应变场出现周期性条纹现象应变云图上出现规则的正负交替条纹像水波纹。原因step_size和strain_window不匹配差分窗口内点数太少位移的整像素量化误差被放大成周期性噪声。解决把strain_window增大到 9 或 11或者把step_size减小到 3让差分窗口覆盖更多点。代价是空间分辨率下降需要根据关心的最小应变尺度权衡。4.4 内存溢出或程序越跑越慢现象跑到几十帧后 MATLAB 卡死或报Out of memory。原因循环里不断累积变量或者parfor并行时每个 worker 都复制了一份大图像。解决在帧循环末尾显式clear临时变量用pack整理内存碎片并行计算时把图像序列改成按需读取不要一次性imread所有帧。如果用了parfor注意params结构体要设成广播变量避免每个 worker 重复拷贝。4.5 应力结果和试验机读数对不上现象DIC 算出的平均应力和拉伸机载荷除以截面积的结果差很多。原因pixel_size标定错误或者 ROI 区域选得不对把非均匀变形区比如夹持段也算进去了。解决用标定板重新标定pixel_size确保 ROI 只覆盖均匀变形的工作段。另外检查本构参数E和nu是不是用的实际材料值不要用默认的钢参数去算铝合金。5. 二次开发与改造把这套代码变成自己的实验工具5.1 换数据源从 TIFF 序列到视频流默认的load_sequence.m读的是 TIFF 序列如果你的数据是 MP4 或相机实时流改这个文件就行。核心是把imread换成VideoReader的readFrame保持输出格式一致灰度、double 或 uint8 都行但整个流程要统一。function frames load_sequence(src, type) if strcmp(type, tiff) files dir(fullfile(src, *.tif)); frames arrayfun((f) imread(fullfile(src, f.name)), files, ... UniformOutput, false); elseif strcmp(type, video) vr VideoReader(src); frames {}; while hasFrame(vr) frames{end1} rgb2gray(readFrame(vr)); %#okAGROW end end end注意rgb2gray只在彩色视频时需要如果是灰度视频直接readFrame即可。frames用 cell 数组存是为了兼容不同尺寸的帧但 DIC 要求所有帧尺寸一致所以实际使用时最好在读取后加一个尺寸校验。5.2 加一个批量处理入口单次实验往往有多组试样手动改路径跑几十次不现实。我一般会在main.m外面再包一层batch_run.m扫描数据目录下的子文件夹每个子文件夹当一组试验自动生成结果目录。function batch_run(root_dir) cases dir(root_dir); cases cases([cases.isdir] ~startsWith({cases.name}, .)); for k 1:numel(cases) case_dir fullfile(root_dir, cases(k).name); fprintf(处理: %s\n, cases(k).name); try main(case_dir); % main 接受数据目录参数 catch ME fprintf(失败: %s - %s\n, cases(k).name, ME.message); end end endtry-catch是必须的一组数据出错不能影响后面的。main函数改成接受目录参数而不是硬编码路径这样批量调用才成立。结果目录按试样名自动创建避免覆盖。5.3 结果验证用合成图做回归测试改完算法后怎么确认没改坏最可靠的办法是生成一张已知位移场的合成图跑一遍看输出和理论值的误差。我习惯在tests/目录下放一个test_synthetic.m每次改完核心算法先跑它。% 生成已知正弦位移场的合成散斑图 [x, y] meshgrid(1:512, 1:512); speckle rand(512); u_true 2 * sin(2*pi*x/256); % 水平位移幅值 2 像素 v_true 1.5 * cos(2*pi*y/256); % 垂直位移幅值 1.5 像素 % 用 interp2 做图像变形 def interp2(speckle, xu_true, yv_true, linear, 0); % 跑 DIC比较输出和 u_true/v_true [u_est, v_est] dic_corr(speckle, def, params); err_u u_est - u_true(1:params.step_size:end, 1:params.step_size:end); fprintf(水平位移 RMS 误差: %.4f 像素\n, rms(err_u(:)));这个测试的价值在于位移场是解析已知的误差可以直接量化。如果 RMS 误差在 0.05 像素以内说明整像素加亚像素流程正常如果突然跳到 0.5 像素多半是某次改动引入了 bug。合成图还能测边界情况比如大位移、大梯度、低对比度比真实实验数据更容易定位问题。5.4 性能优化从分钟级到秒级原始双重循环在 2048×2048 图像、step 5 的情况下单帧可能要跑几分钟。几个实用的加速手段一是把内层循环向量化用im2col把子区提取变成矩阵操作二是用parfor替换外层for按计算点并行三是把normxcorr2换成频域 FFT 实现大子区时快很多。我一般先做parfor改动最小、收益最明显通常能提速 4~8 倍具体看核数。从那以后我每次改完core/里的任何文件都强制先跑一遍合成图回归测试确认 RMS 误差没退化再上真实数据。这个习惯帮我挡掉了至少三次“看起来能跑但结果全错”的改动。希望这套拆解能帮你少走点弯路把 DIC 真正用起来。本文还有配套的精品资源点击获取