
简介基于MATLAB的瞳孔定位资源以阈值分割和灰度分布特性为技术主线适合生物识别、医学诊断、人机交互等方向的图像处理学习者与开发者用于解决复杂背景中瞳孔区域的自动提取与边缘定位问题。压缩包共5个文件包括4个M脚本和1张JPG示例图像脚本覆盖图像读取、直方图分析、阈值二值化、形态学优化及边界提取等环节示例图像可验证完整流程。已有447人浏览学习。通过源码可掌握imbinarize、histcounts、bwboundaries等函数的实际组合用法理解依据灰度直方图波峰波谷选取阈值的原理并学习应对光照变化、轮廓不完整等情况的形态学修正与迭代策略从灰度分布分析到瞳孔中心与边缘确定均提供了可运行的参考实现。1. 阈值分割与瞳孔定位为什么灰度分布比边缘检测更稳在虹膜识别和眼动追踪里瞳孔定位最常见的失败场景是眼皮遮挡产生乱边缘、镜片反光造成高亮斑块、睫毛阴影形成伪轮廓。这时候单纯用 canny 或 sobel 找边缘往往会把眼睑和反光点当成瞳孔边界。反过来看灰度分布瞳孔区域在近红外或普通摄像头上都稳定地是眼球里最暗的连通块虹膜和巩膜灰度明显更高。利用这个先验先用阈值把图像切分成暗区和亮区再靠连通域筛选把瞳孔捞出来比边缘检测更抗干扰。pupil-localization 压缩包里包含 MAINCODE.m 主程序、binaryzation.m 二值化、drawcircle.m 画圆、compiris.m 灰度分布比较以及示例图像 04L.jpg恰好像一条从直方图分析到瞳孔圆拟合的链路适合 matlab 图像处理课程设计或入门级瞳孔定位落地。2. 灰度直方图分析与阈值选择从双峰到 Otsu2.1 瞳孔低灰度先验为什么直方图能直接反映瞳孔位置瞳孔在光学上是一个近似黑体的开孔进入眼球的光线大部分被吸收所以成像后瞳孔区域的灰度值通常落在整个图像的最低 10%20% 区间。虹膜和巩膜的反射率远高于它因此在灰度直方图上会形成一个明显的主峰巩膜和一个次峰或不规则缓坡虹膜瞳孔则以小包或者长尾形式出现在低灰度侧。需要说明的是普通可见光图像里瞳孔的灰度峰值有时很低只有几个到几十个灰阶如果使用近红外图像瞳孔反而可能更亮这时判断逻辑要反向。本项目的 04L.jpg 是典型可见光人眼图按暗前景处理即可。用 imhist 直接看直方图是一种最省事的读图方式但它只输出 counts没有把灰阶与像素数一一对应成变量。后面还要自动找谷底我一般会用 histcounts 拿边界和计数再配合 movmean 做平滑避免把噪声尖峰误判成模板峰。2.2 用 histcounts 计算灰度分布并定位谷底% 读取示例图像并转灰度 I imread(04L.jpg); if size(I, 3) 3 I rgb2gray(I); end % 统计 0~255 每个灰阶的像素数 [counts, edges] histcounts(I, BinLimits, [0 255], BinWidth, 1); grayLevel (edges(1:end-1) edges(2:end)) / 2; % 平滑计数防止伪峰 smoothed movmean(counts, 5); % 找主峰和次峰取峰间谷底作为分割阈值 [pkVals, pkLocs] findpeaks(smoothed, MinPeakHeight, 50, MinPeakDistance, 20); if length(pkLocs) 2 [~, betweenIdx] min(smoothed(pkLocs(1):pkLocs(2))); thresholdVal grayLevel(pkLocs(1) betweenIdx - 1); else thresholdVal graythresh(I) * 255; endhistcounts 的BinLimits指定统计范围BinWidth1让每个灰阶一个 bingrayLevel取每个区间的中点这样之后得到阈值后可以直接做像素比较。movmean是对计数值做窗口为 5 的移动平均把相机噪声造成的小锯齿抹平。findpeaks用MinPeakHeight过滤掉零星高峰用MinPeakDistance20保证两个峰至少相隔 20 个灰阶避免把同一峰的抖动识别成两个峰。峰之间最小计数值对应的灰阶就是视觉上的谷底阈值取在这里比取在峰腰更稳。findpeaks需要 Signal Processing Toolbox没有的话可以用循环找峰替代后面 MAINCODE.m 里会给出一个不依赖该工具箱的版本。2.3 从谷底到切割二值化前的灰度分布验证在真正做二值化之前我建议先把谷底阈值放到直方图里画出来确认它落在瞳孔包和虹膜峰之间的凹陷处。如果阈值选得太靠左瞳孔会被切掉一部分二值化后瞳孔面积偏小太靠右会把虹膜深色区域也划进来产生大面积粘连。这里给出一个可视化片段figure; histogram(I, BinLimits, [0 255], BinWidth, 1, FaceAlpha, 0.3); xline(thresholdVal, r-, threshold);histogram 用同样的 bin 设置保证与前面统计一致xline画出阈值位置。这一步不是核心计算但在调参时能省很多时间尤其当 04L.jpg 之外的图像光照不一致时。肉眼看到阈值左侧的像素基本都集中在黑色瞳孔区域右侧保留虹膜纹理说明谷底位置找得对。2.4 固定阈值、Otsu 与 isodata 的取舍谷底法本质是固定阈值。固定阈值对同一相机、同一光照环境下的批量图像足够用但换一个场景就需要重新标定。Otsu 通过最大化类间方差自动选择一个阈值适合直方图近似双峰的场景缺点是瞳孔区域像素占比很小或光照不均匀时Otsu 会把虹膜暗部并入瞳孔。isodata 迭代法以当前阈值为界计算两侧均值再取均值的中点作为新阈值收敛后得到的结果更偏向把暗部中心保住。方法依赖本场景风险适用谷底法直方图双峰明显只对固定光照有效课程设计、批量固定环境Otsu类间方差最大瞳孔占比小时易过分割快速尝试验证isodata迭代收敛收敛慢需设定初值光照有缓慢漂移Otsu 在 MATLAB 里就是一行graythresh(I)它返回归一化阈值。isodata 可以用自写循环实现这里不展开。实际工程里我常用两个阈值做带通thrLow I thrHigh这样既排除比瞳孔还暗的镜架深色区域也排除虹膜高光比单阈值多一层保护。这个技巧会在后面的 MAINCODE.m 流程里体现。3. MAINCODE.m 与 binaryzation.m二值化与连通域筛选3.1 binaryzation.m 函数设计与参数边界binaryzation.m 从名字看就是封装二值化步骤。这里给出一个与项目语义一致的常见实现输入灰度图输出逻辑矩阵true 表示瞳孔候选像素。为了让不同图像都能跑函数内部把阈值写成可配置参数缺省时用 Otsu 兜底function bw binaryzation(I, threshold) % I: 灰度图像 % threshold: 可选标量灰阶阈值 if nargin 2 threshold graythresh(I) * 255; end bw I threshold; end这里用I threshold而不是imbinarize(I, threshold/255)是为了避免 imbinarize 在部分旧版本里对 uint8 输入做额外规范化。需要注意 threshold 的单位是灰阶 0~255而 graythresh 返回的是 0~1 归一化值所以乘以 255。如果原项目里的 binaryzation.m 用自己的统计逻辑比如结合局部均值效果差异主要在光照剧烈变化时体现参数含义不变。3.2 形态学开闭运算清洗噪声并保持瞳孔面积二值化后瞳孔区域内可能有睫毛造成的孔洞区域外有反光形成的小白斑。形态学闭运算能先填充孔洞开运算再删除细微突起。结构元素用圆盘因为瞳孔本身接近圆形圆盘结构元素不会引入方向偏差。se strel(disk, 5); % 半径 5 像素对应 04L.jpg 瞳孔半径约 40 像素 bw imclose(bw, se); % 先闭运算填充孔洞 bw imopen(bw, se); % 再开运算去掉边缘毛刺strel(disk, 5)里的 5 是结构元素半径需要根据图像分辨率同比缩放。我的经验是取瞳孔期望半径的 1/10 左右太小挡不住睫毛太大会把瞳孔边缘磨圆。下面的表格给出调整方向结构元素半径效果副作用过小孔洞补不上边界保留细碎适中闭运算填充孔洞开运算平滑边界基本无副作用过大边界外扩真实瞳孔面积被侵蚀MATLAB 中 r2017a 之后 strel 使用也保持兼容matlab r2023b 及更早版本都能直接运行这段代码。3.3 在 MAINCODE.m 中串联直方图、二值化与大连通域提取MAINCODE.m 是主脚本。常规做法是读图、转灰度、计算直方图谷底阈值、调用 binaryzation然后通过 regionprops 把最大连通域作为瞳孔候选。这个流程把上一章的直方图分析落到实处。I imread(04L.jpg); if size(I, 3) 3 gray rgb2gray(I); else gray I; end % 步骤 1: 确定灰度阈值这里用 Otsu 初值可替换为第 2 章谷底法结果 thresholdVal graythresh(gray) * 255; % 步骤 2: 调用 binaryzation 得到暗区 bwPupil binaryzation(gray, thresholdVal); % 步骤 3: 形态学清洗 bwPupil imclose(bwPupil, strel(disk, 5)); bwPupil imopen(bwPupil, strel(disk, 5)); % 步骤 4: 提取最大连通域 stats regionprops(bwPupil, Area, Centroid, BoundingBox); [~, maxIdx] max([stats.Area]); cx stats(maxIdx).Centroid(1); cy stats(maxIdx).Centroid(2); pupilBox stats(maxIdx).BoundingBox;regionprops 返回的是一个结构数组每个连通域一个元素。[stats.Area]把所有面积拼成向量max 找到最大索引然后取质心作为瞳孔中心。这里有两个隐藏问题其一如果瞳孔没有与镜框暗区完全分离最大连通域可能是镜框需要用面积范围和圆形度过滤其二质心是像素质心不是圆几何中心瞳孔边缘缺失时质心会偏移。常见做法是再结合Circularity属性过滤circularities [stats.Circularity]; candidates find([stats.Area] 1000 circularities 0.6);这一步我一般放在 max 之前避免选中长条形睫毛区域。MAINCODE.m 里如果没有这行你需要根据自己的图像决定是否加上。注意regionprops 的 Centroid 返回的是 [x y] 格式x 是列方向y 是行方向。如果你习惯用矩阵下标 [row col]在绘图前要交换顺序。4. drawcircle.m 与 compiris.m边界拟合与灰度分布比较4.1 drawcircle.m 的标准画法drawcircle.m 在项目里的作用是把拟合出的瞳孔圆叠加到原图上方便人眼复核。它通常接收圆心坐标、半径和坐标轴句柄用参数方程画 100 个点的多段线。function h drawcircle(ax, cx, cy, r) theta linspace(0, 2*pi, 100); x cx r * cos(theta); y cy r * sin(theta); h plot(ax, x, y, r-, LineWidth, 2); endlinspace(0, 2*pi, 100) 生成 100 个采样角采样越密圆越平滑100 已经足够显示。调用时只需drawcircle(gca, cx, cy, r)。注意这里画的是叠加线不是掩膜所以它不影响后续像素操作。半径 r 可以从 regionprops 得到的 Area 计算等效半径sqrt(area/pi)也可以用 BoundingBox 的宽高平均再除以 2。对于非正圆瞳孔drawcircle 会丢掉椭圆倾斜信息但如果只需要中心点这个近似足够。4.2 compiris.m 的灰度分布比较逻辑compiris 这个名字可能是 compare iris 的缩写也可以理解为 compare inner/outer region intensity。它在项目里大概率是验证分割质量瞳孔内部灰度应该显著低于四周虹膜区域。用极坐标或环形掩膜实现。% 生成瞳孔内部掩膜 [rows, cols] size(gray); [xx, yy] meshgrid(1:cols, 1:rows); maskIn (xx - cx).^2 (yy - cy).^2 r^2; % 生成瞳孔外环掩膜: 半径 r 到 1.4r maskRing (xx - cx).^2 (yy - cy).^2 r^2 ... (xx - cx).^2 (yy - cy).^2 (1.4*r)^2; meanIn mean(gray(maskIn)); meanRing mean(gray(maskRing)); contrastRatio meanRing / (meanIn eps);maskIn 是逻辑矩阵元素为 true 的位置属于瞳孔内gray(maskIn)用逻辑索引把所有满足条件的像素取出来。meanRing 计算虹膜内环的灰度均值。contrastRatio 大于 1 说明瞳孔确实比周围暗。eps加在分母上防止瞳孔内部全黑导致除零。这个比值可以直接作为定位置信度低于 1.1 时说明分割把虹膜或阴影误判成了瞳孔需要调阈值。场景contrastRatio 典型值判断正常瞳孔 1.3定位可信瞳孔反光1.0~1.2边缘可能偏移需回看二值图严重眼睑遮挡 1.0分割失败改用边缘约束4.3 光照补偿与鲁棒性处理瞳孔定位真正难处理的是光照不均。这时候固定阈值在图像一侧失效compiris 的比值也会变低。常规手段是先做直方图均衡化再用 adaptive threshold最后做形态学。MATLAB 里可以用adapthisteq做对比度受限的自适应直方图均衡化它比histeq更适合局部光照变化。这一步我会放在灰度化之后、直方图统计之前避免把反射高光放大成斑块。要注意如果项目原图是红外图瞳孔可能比虹膜亮那么 binaryzation 里的比较符号要反转compiris 的分子分母也要交换。代码里应把gray threshold写成可配置的darkPupil标志。这部分不要写死否则换摄像头就废了。另外04L.jpg 如果带 EXIF 方向的旋转信息直接 imread 后可能出现行列倒置建议先用imfinfo检查 Orientation 字段必要时用imrotate修正再继续。5. 用 04L.jpg 实测参数调整与验证方法5.1 跑通 MAINCODE.m 的环境与文件路径04L.jpg 放在 pupil-localization 根目录确保 MAINCODE.m 的工作路径指向该目录。如果输入图像是 uint8 RGBrgb2gray 后直接参与 histcounts 不会报错但有些脚本会把图像归一化到 0~1 再传参这样阈值也会变到 0~1需要保持一致。建议在 MAINCODE.m 开头加一句assert(isa(gray,uint8), 请先转换到 uint8 灰度图)避免单位混乱。5.2 中间结果可视化与参数排查运行时代码后面加一个 subplot 窗口把原图、直方图、二值图、拟合圆放在一起看figure; subplot(2,2,1); imshow(gray); title(gray); subplot(2,2,2); histogram(gray, BinLimits, [0 255], BinWidth, 1); hold on; xline(thresholdVal,r-); subplot(2,2,3); imshow(bwPupil); subplot(2,2,4); imshow(gray); hold on; drawcircle(gca, cx, cy, r);如果二值图里瞳孔区域有孔把闭运算的 disk 半径加大如果瞳孔与背景粘连增大阈值或改用双阈值带通。常见问题如下现象原因调整二值图瞳孔区域偏小阈值太低增大 thresholdVal瞳孔与虹膜粘连阈值太高减小 thresholdVal 或改用带通质心偏向睫毛未过滤低圆形度连通域使用 Circularity 过滤拟合圆明显偏离瞳孔被反光剖开先闭运算再开运算5.3 用质心偏差定量验证定位精度最后给一个简单有效的验证技巧手动在图像上点出瞳孔中心然后计算算法质心与手动中心之间的欧氏距离。这个指标比肉眼看圆更客观也容易写进实验报告。% 手动标记瞳孔中心比如用 ginput 或直接选定坐标 manualCenter [x_m, y_m]; % 从图上观察得到 pupilCenter [cx, cy]; centerErr sqrt(sum((pupilCenter - manualCenter).^2));centerErr 小于 3 像素通常说明定位准确3~10 像素说明有轻微偏差可接受大于 10 像素就需要检查阈值和形态学参数。如果你想更严格可以把半径也纳入验证手动点瞳孔左右和上下的四个边缘点用最小二乘拟合圆心与半径再与 regionprops 结果比较这样能同时评估边界误差。本文还有配套的精品资源点击获取