简介这是一套基于MATLAB实现的数字图像相关DIC主程序资源面向材料科学、工程力学与生物医学领域的学习者与研究人员用于通过连续图像序列非接触式测量物体表面的位移、应变与应力。资源包共161个文件包含82个m脚本、43张jpg与19张png图像、10个html说明页、5个fig图形文件及少量txt文档压缩包约6.01MB代码与示例图像、界面文件配套齐全。内容覆盖图像预处理、角点检测与特征配对、互相关匹配、位移场计算、应变应力求解及结果可视化等完整流程脚本注释与文件结构便于逐步调试运行。目前已有4164人学习下载适合初学者理解DIC算法原理也方便进阶用户按实验条件修改优化代码快速搭建自己的光学测量分析流程。1. 从一叠试件照片到全场应变云图DIC 主程序到底替你省掉了什么手头有一组拉伸试验的序列照片试件表面喷了黑白散斑相机固定、光源稳定几十到几百张图按时间顺序排好。你想知道的不是“试件断没断”而是每一帧之间材料表面每个点往哪个方向动了多少、局部应变集中在哪里、应力是怎么随着位移场演化的。拿卡尺量、拿引伸计量只能得到几个离散点的数据而数字图像相关DICDigital Image Correlation要做的是把整幅图像变成一个密集的位移场再从这个位移场里算出应变场最后结合材料本构关系推出应力分布。这套 DIC MATLAB 主程序解决的就是“从图像序列到位移、应变、应力”这条链路。它不负责采集图像也不负责喷散斑它负责的是最核心也最容易翻车的一步把序列图像读进来划好感兴趣区域ROI做子区匹配算出整场的 U、V 位移再通过差分或最小二乘拟合得到应变最后按你给定的弹性模量和泊松比换算成应力。适合谁用做材料力学实验的研究生、做结构变形监测的工程师、以及想用 MATLAB 快速验证 DIC 算法而不是从零写匹配代码的人。如果你手里已经有序列图缺的是一个能跑通、能调参、能出云图的主程序那这份资源就是冲着你来的。2. 主程序怎么读图、划 ROI、做子区匹配三个环节的参数逻辑2.1 图像序列的读取与预处理别让第一帧就埋雷DIC 主程序的第一步不是算而是把图读对。常见做法是用imread配合dir批量读取但这里有个容易被忽略的点文件名排序。如果相机导出的是img1.tif、img2.tif……img10.tif直接按字符串排序会变成img1、img10、img2序列顺序全乱后面算出来的位移就是跳变的。我一般会在读取后强制按文件名中的数字重新排一次。% 读取指定文件夹下的图像序列并按文件名中的数字排序 imgFolder D:\DIC_test\specimen_01; % 图像所在文件夹 imgFiles dir(fullfile(imgFolder, *.tif)); % 按需改成 *.png / *.jpg % 提取文件名中的数字用于排序避免 img10 排在 img2 前面 numIdx zeros(length(imgFiles), 1); for k 1:length(imgFiles) token regexp(imgFiles(k).name, \d, match); numIdx(k) str2double(token{end}); end [~, sortOrder] sort(numIdx); imgFiles imgFiles(sortOrder); % 读取第一帧作为参考图并转为灰度 refImg imread(fullfile(imgFolder, imgFiles(1).name)); if size(refImg, 3) 3 refImg rgb2gray(refImg); end refImg im2double(refImg);这段代码的逻辑说明dir拿到的是文件结构体数组regexp提取文件名里最后一组数字作为排序依据sort返回的sortOrder才是真正的时间顺序。参数上imgFolder换成你自己的路径文件扩展名按实际格式改。转灰度是因为 DIC 匹配通常在灰度图上做彩色图反而增加计算量且不提升精度。im2double把像素值归一化到 0 到 1后续做归一化互相关时数值更稳定。预处理里还有一个常被跳过但影响很大的步骤ROI 的划定。主程序一般会提供一个掩膜或者矩形区域让你只对试件表面做计算避免背景区域参与匹配产生伪相关。我一般会先用第一帧画一个矩形 ROI把试件边缘裁掉几个像素因为边缘子区容易超出图像边界导致匹配失败。2.2 子区匹配的核心参数子区大小、步长、搜索窗口DIC 的位移计算本质是在参考图上取一个子区subset在变形图上滑动一个搜索窗口找归一化互相关NCC或平方差SSD最大的位置那个位置偏移就是位移。主程序里通常暴露三个关键参数子区大小subsetSize、计算步长stepSize、搜索半径searchRadius。子区大小决定了空间分辨率和噪声抑制的平衡。子区太小比如 11×11 像素对散斑质量要求极高噪声一大匹配就飘子区太大比如 61×61位移场会被过度平滑应变集中区域算不出来。常见做法是取 21 到 41 之间散斑颗粒直径大约占 3 到 5 个像素时子区取 31×31 是比较稳的起点。步长决定输出位移场的密度。步长等于子区大小时输出的是不重叠的网格点步长小于子区大小输出的是重叠的密集场。做应变分析时我一般会把步长设成子区大小的 1/4 到 1/2比如子区 31、步长 8这样应变计算有足够的数据点做局部拟合。搜索半径要覆盖你预估的最大位移。如果两帧之间试件移动了 20 个像素搜索半径至少给 25 到 30。给太小匹配直接跑到窗口边缘位移被截断给太大计算量上去而且可能匹配到错误位置。一个实用技巧是先用大搜索半径跑前两帧看最大位移量级再回头收紧。% DIC 主程序核心参数设置 subsetSize 31; % 子区边长奇数常用 21~41 stepSize 8; % 计算步长越小场越密计算量越大 searchRadius 30; % 搜索半径需大于预估最大位移 nccThreshold 0.85; % 相关系数阈值低于此值视为无效点 % 在参考图上按步长生成计算网格点 [imgH, imgW] size(refImg); halfSub floor(subsetSize / 2); gridX halfSub 1 : stepSize : imgW - halfSub; gridY halfSub 1 : stepSize : imgH - halfSub; [gridX, gridY] meshgrid(gridX, gridY); pointsRef [gridX(:), gridY(:)]; % 每个点的参考坐标参数说明nccThreshold是过滤无效匹配的阀门低于这个相关系数的点直接标记为 NaN不参与后续应变计算。这个阈值不要设太低0.7 以下基本是噪声在主导也不要设太高0.95 以上在散斑质量一般时会把大量有效点误杀。0.85 到 0.9 是常见区间。2.3 从位移场到应变场差分、拟合与应力换算拿到每个网格点的 U、V 位移后应变不是简单地对相邻点做差分就完事。位移场里含有匹配噪声直接差分会把噪声放大成应变的剧烈波动。主程序里通常提供两种方式一种是局部最小二乘拟合在每一个点周围取一个小窗口拟合一个平面平面的斜率就是应变另一种是全局平滑后再差分。我一般会优先用局部最小二乘窗口大小取 5×5 或 7×7 个计算点。窗口太小噪声压不住窗口太大应变集中区域被抹平。这个窗口和子区大小是独立的子区管匹配拟合窗口管应变平滑。应力换算相对直接在弹性阶段用胡克定律平面应力状态下 σx E/(1-ν²) (εx ν εy)σy 同理τxy G γxy。主程序里会让你输入弹性模量 E 和泊松比 ν然后逐点算出应力分量。这里有个前提这套换算只适用于线弹性、小变形。如果你的试件已经进入塑性或者大变形应力结果只能当参考不能当定量结论。% 局部最小二乘拟合计算应变场以 ex 为例 strainWindow 5; % 拟合窗口半径单位是计算点 E 210e3; % 弹性模量单位 MPa nu 0.3; % 泊松比 h stepSize; % 计算点间距对应物理尺寸需按标定换算 % 假设 U 为位移场矩阵尺寸与 pointsRef 网格一致 % 对每个有效点取其邻域做平面拟合 [ny, nx] size(U); ex nan(ny, nx); ey nan(ny, nx); gxy nan(ny, nx); for i 1 strainWindow : ny - strainWindow for j 1 strainWindow : nx - strainWindow % 取邻域位移 uPatch U(i-strainWindow:istrainWindow, ... j-strainWindow:jstrainWindow); vPatch V(i-strainWindow:istrainWindow, ... j-strainWindow:jstrainWindow); if any(isnan(uPatch(:))) || any(isnan(vPatch(:))) continue; % 含无效点的邻域跳过 end % 平面拟合求偏导此处用简单差分近似示意 [du_dx, ~] gradient(uPatch, h); [~, dv_dy] gradient(vPatch, h); [du_dy, ~] gradient(uPatch, h); [dv_dx, ~] gradient(vPatch, h); ex(i,j) mean(du_dx(:)); ey(i,j) mean(dv_dy(:)); gxy(i,j) mean(du_dy(:) dv_dx(:)); end end % 平面应力状态下的应力换算 sigmaX E / (1 - nu^2) * (ex nu * ey); sigmaY E / (1 - nu^2) * (ey nu * ex); tauXY E / (2 * (1 nu)) * gxy;逻辑说明gradient在这里做的是数值偏导h是计算点之间的物理距离需要你根据标定板或者试件已知尺寸换算。如果h用的是像素单位那算出来的应变是“每像素”的应变不是真实应变这一步是很多人翻车的地方。应力换算的三个公式对应平面应力线弹性E和nu按材料实际值填。算完之后无效点仍然是 NaN画云图时会自动留白。3. 避坑与排查DIC 主程序跑不出结果的五个血泪经验3.1 位移场出现大面积 NaN相关系数普遍偏低现象跑完主程序位移云图上一大片空白检查nccThreshold发现大量点低于阈值。原因通常不是算法错而是散斑质量或者光照变化。散斑颗粒太大、太小、对比度不够或者序列图像之间光源闪烁导致灰度整体漂移都会让 NCC 峰值不明显。解决先单独看第一帧和第十帧的直方图如果均值差超过 10 个灰度级说明光照在变需要做灰度归一化或者重新采集。散斑方面颗粒直径占 3 到 5 个像素是经验值喷得太密或太稀都要重做。3.2 位移场整体偏移一个固定值像是“零点漂移”现象刚体平移试验里试件没动但算出来的 U 场有一个非零常数。原因多半是参考子区和变形子区的坐标对应关系搞错了比如网格点生成时多加了半个子区或者图像读取时第一帧不是真正的参考帧。解决用已知刚体平移量做标定比如手动把第二帧整体移 5 个像素看算出来是不是 5。如果不是检查pointsRef的生成逻辑和imread的帧顺序。3.3 应变云图噪声极大像雪花屏现象位移场看起来还算平滑但应变场全是高频抖动。原因是位移场里的匹配噪声被差分放大了。解决两条路一是增大子区尺寸从 31 加到 41 或 51用空间换噪声二是增大应变拟合窗口从 5 加到 7 或 9。如果还不行说明散斑质量本身不够匹配噪声底子太大再怎么平滑也是徒劳。3.4 应力结果量级明显不对差几个数量级现象算出来的应力是几百 MPa 或者零点几 MPa和材料实际强度对不上。原因几乎都是单位问题。位移场如果是像素单位应变就是无量纲的“像素应变”但如果你在应力公式里用的E是 MPa而h没有换算成毫米那应力就全错了。解决先确认标定比例比如 1 像素对应多少毫米把h换成物理单位再检查E的单位是 MPa 还是 Pa。我一般会在主程序开头强制写一行单位注释跑之前先读一遍。3.5 计算速度慢到无法接受一帧要跑几分钟现象网格点几千个每个点都在全图搜索MATLAB 循环扛不住。原因是搜索策略太暴力没有用频域加速或者预测搜索。解决常见做法是用normxcorr2做模板匹配它底层是频域实现比手写双重循环快很多另一个技巧是用上一帧的位移作为下一帧的搜索中心预测把搜索半径缩小到实际位移的 1.5 倍左右。如果还是慢把步长从 8 调到 16计算点数量直接降到四分之一。4. 进阶用法批量处理、位移场后处理与结果验证4.1 批量跑完整个序列自动保存每帧的位移和应变单帧跑通之后下一步是把整个序列串起来。主程序一般会提供一个循环框架但这里有个细节参考帧的选择。常见做法有两种一种是固定第一帧为参考帧后续每一帧都和第一帧匹配适合小变形累积另一种是相邻帧匹配适合大变形但需要把增量位移累加。我一般会先用固定参考帧跑一遍如果后面帧的相关系数掉得厉害再切到增量模式。% 批量处理整个序列固定第一帧为参考 nFrames length(imgFiles); U_all cell(nFrames, 1); V_all cell(nFrames, 1); for f 1:nFrames curImg imread(fullfile(imgFolder, imgFiles(f).name)); if size(curImg, 3) 3 curImg rgb2gray(curImg); end curImg im2double(curImg); % 调用主程序的核心匹配函数返回 U、V 场 [U, V, nccMap] dic_match(refImg, curImg, ... pointsRef, subsetSize, searchRadius, nccThreshold); U_all{f} U; V_all{f} V; fprintf(Frame %d / %d done, valid points: %d\n, ... f, nFrames, sum(~isnan(U(:)))); end % 保存结果避免重复计算 save(dic_results.mat, U_all, V_all, pointsRef, ... subsetSize, stepSize, searchRadius);逻辑说明dic_match是主程序里的核心函数输入参考图、当前图、参考点坐标和三个参数输出位移场和相关系数图。fprintf打印每帧的有效点数方便你实时判断有没有帧跑崩。最后用save把结果存成 MAT 文件后面调参重算时不用再跑一遍匹配。4.2 用刚体平移和已知应变做验证别直接信云图DIC 结果看起来漂亮不代表对。我一般会做两个验证第一刚体平移验证把一张图人为平移整数像素跑 DIC看位移场是不是均匀等于平移量误差应该在 0.01 到 0.05 像素量级第二已知应变验证用模拟散斑图施加一个已知的线性位移梯度看算出来的应变和设定值差多少。这两个验证跑通再去看真实试验的云图心里才有底。验证项目输入期望输出可接受误差刚体平移图像整体移 5 pxU 场均匀等于 5±0.05 px均匀应变线性位移梯度应变场常数±50 με无效点过滤遮挡区域NaN 标记不参与应变4.3 位移场后处理平滑、插值和应力集中识别原始位移场总带噪声直接算应力会放大。我一般会在应变计算前做一次中值滤波或者高斯平滑窗口取 3 到 5 个计算点。插值方面如果网格点不够密可以用griddata插到更细的网格再画云图但注意插值不产生新信息只是让图好看。应力集中识别靠的是应变梯度我通常会把 von Mises 应力或者最大剪应变单独画出来看集中区域是不是和试件几何缺口或者裂纹尖端对得上。对不上先回头查 ROI 和标定别急着下结论。从那以后我每次跑 DIC 之前都强制先用刚体平移图走一遍全流程确认位移场没有系统性偏移、应变场没有异常噪声才敢把真实试验的序列图丢进去。这个习惯帮我省掉了至少三次把错误结果写进报告的尴尬。希望帮到你。本文还有配套的精品资源点击获取