基于灰度直方图与阈值分割的瞳孔定位MATLAB实现
2026/9/14 15:02:58 网站建设 项目流程

简介:基于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. 灰度直方图分析与阈值选择:从双峰到 Otsu

2.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; end

histcounts 的BinLimits指定统计范围,BinWidth=1让每个灰阶一个 bin,grayLevel取每个区间的中点,这样之后得到阈值后可以直接做像素比较。movmean是对计数值做窗口为 5 的移动平均,把相机噪声造成的小锯齿抹平。findpeaksMinPeakHeight过滤掉零星高峰,用MinPeakDistance=20保证两个峰至少相隔 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); end

linspace(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 RGB,rgb2gray 后直接参与 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 结果比较,这样能同时评估边界误差。

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

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

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

立即咨询