基于随机Hough变换圆检测算法的MRI质量控制空间线性度自动检测程序设计

林迪逵，吴剑威，李威，付丽媛，熊晖，陈自谦

解放军福州总医院 a. 医学工程科/福建省医学装备管理质量控制中心；b. 医学影像中心，福建 福州 350025

引言

MRI扫描结果的精确与否是评判MRI技术临床应用成果的基础标准，我们希望MRI图像能够对患者的解剖、病理和生理功能作出如实呈现，但是问题在于无论对于临床研究或是诊断，实现高质量图像的获取都是一项艰巨的任务。基于上述临床应用背景，MRI研究与临床工作者提出了质量控制（Quality Control，QC）概念，旨在对磁共振成像质量实现监测，并通过一系列的技术程序保证MRI图像质量的持续提高。MRI质量控制程序围绕图像质量参数的检测进行展开，通过成像质量来客观反映成像设备的运行状态，实现对成像系统性能状态的评价。质量控制程结果帮助临床工程师对设备性能的改变作出反应，保证MRI成像系统在安全、稳定的状态下运行[1]。

MRI QC程序对于空间线性度的检测目的是为验证MRI系统是否以反映研究对象真实尺寸的缩放方式进行成像。良好的空间线性度保证MRI图像对成像目标距离进行真实再现，实现对组织几何关系的精确反映，有助于临床医师进行准确的影像诊断。影响磁共振成像空间线性度的因素有很多，不合适的梯度缩放校正因子、不均匀磁场等成像环境都能导致MRI图像空间线性度的下降[2-3]。

确定空间线性度的传统方法依靠人工检测，其存在以下缺陷：① 涉及多组重复测距，过程繁琐耗时；② 检测精度受限于QC检测人员感觉器官的生理变化；③ 若检测人员发生更替，QC的连续性意义也将大打折扣。为解决上述问题，本文拟设计空间线性度检测算法，取代传统人工检测，实现空间线性度的自动测量，保证参数检测结果如实反映成像设备相关性能状态。

1 材料与方法

1.1 Magphan SMR 170体模的空间线性度检测方法

Magphan SMR 170体模由美国体模实验室研制，该体模在横断位上设计5个测试层面，层面上设置的多组结构模块允许MRI QC程序对基本质量控制参数进行检测，达到校验成像系统的目的。

QC程序对空间线性度的检测在体模第四层横断位图像上进行（图1）。该过程主要针对测试层面上若干圆孔间的间距进行测量，定义空间线性度l为：

其中LR为体模中实际圆孔间距，三组孔间实际距离分别为80、40、20 mm；LM为人工检测距离[4-5]。

图1 基于Magphan SMR 170体模的空间线性度检测示例

从上述检测过程可以看出，若要实现空间线性度的自动检测，关键在于对图像中若干目标模块（即圆孔）的识别，同时圆孔定位的精确性决定了线性度测量结果的准确性。

1.2 基于随机Hough变换的圆检测

在数字图像的自动检测过程中，对目标形状的检测要求常常作为子问题出现，例如实现对直线、圆或椭圆的识别。在大多数情况下，在预处理阶段通过一些边缘检测程序的实施，可以获得图像空间中期望曲线上的图像点或图像像素。但是由于图像数据质量较差或者边缘检测程序存在缺陷等原因，可能存在期望曲线上的点或像素缺失等情况，导致上述边缘噪声点与理想的目标形状将存在空间偏差。所以，为达到目标形状的检测要求，将提取到的边缘特征聚类成一系列合适的目标形状就显得尤为重要，霍夫变换（Hough Transform）即为解决上述问题而提出的。Hough变换针对参数化的图像对象，通过设置累加器执行“投票”过程，对峰值点作出判断，从而实现将一系列边缘点列为候选对象的目标。

基于标准Hough变换（Standard Hough Transform，SHT）的圆检测算法由于引进了高维参数空间，在某些图像实时处理的应用环境中，变换效果并不理想。为解决SHT算法缺陷，随机Hough变换（Randomized Hough Transform，RHT）应运而生。RHT的特点在于采用了随机抽样机制，实现图像空间到参数空间的收敛映射，即随机选取图像空间上不共线的3个点，将上述三点映射至参数空间的一个点。相对于SHT的“一到多”映射，RHT提出的“多到一”映射无论对于内存的占用，还是对于时间的消耗，都在很大程度上节省了开销[6-7]。RHT圆检测原理如下。

随机选取待识别圆上不共线的三点（x1，y1）（x2，y2）（x3,y3），分别代入圆方程可得方程组：

以上方程组能确定一个圆心为（a0，b0），半径为r0的圆o0。下面的任务在于确定该待识别圆是否为真实的圆。取第四个圆上点（x4， y4），带入圆心为（a0，b0）的圆方程，得到圆半径1r′。将r0、1r′代入下式：

设 ε0为预设误差值。若 ε1＜ε0，认为点（x4， y4）在圆 o0上，o0确定为候选圆。将待识别圆上的其余点（xi， yi）下的εi值，若εi＜ε0，累加值k加1，当k达到预设阈值kthreshold时，则认为o0为真实的圆。

上述RHT随机采样机制的一个特点在于采样的无目标性，虽然相较于SHT它在算法结构上作出了有效的改进，但是无法避免无效单元的引进（图2）。图2中，o1、o2为待识别的目标圆；d1、d2、d3为三个随机抽样点，其位于不同目标圆上且不共线。根据RHT原理，认为上述3个点可以确定候选圆o3，其圆参数为p，之后将进行搜索参数单元集p等一系列程序来确定圆o3是否为真实圆。由于圆o3为虚假圆，所以后续对于圆o3的检测程序均为无效累积过程。

图2 RHT无效累积示例

假设有N个待识别圆，q为每个待识别圆上所选取点的数量，n为非圆上点。则随机选择的3点落在同一个待识别圆上的概率为：

如果剔除非圆点n，上式简化为：

本文空间线性度待检测图像上待识别的圆孔一共有13个，如果不考虑体模图像中非目标圆上的点，由式(7)可知，随机选择的3点落在同一个待识别圆孔上的概率为1/169，无效累积概率为168/169，如果将非圆边缘点考虑在内，无效累积现象将更为严重。针对本文圆检测的应用环境，认为可以对上述RHT算法作出改进，重点解决无效累积问题，进一步提高圆检测效率。

1.3 改进的RHT圆检测算法

1.3.1 利用梯度方向信息减少无效累积

o1、o2为待识别的目标圆；d1、d2、d3为三个随机抽样点，其位于不同目标圆上且不共线（图3）。由于d1与d2由于位于圆o1上，故d1与d2梯度方向所在的直线与必经过圆心o1，且d1与d2连线的中垂线也必经过圆心o1。反观d1与d3，虽直线与将相交于点d4，但d1与d3连线的中垂线将不再经过点d4，由此可以认定d1与d3不在同一个圆。基于上述判断，后续算法将不再进行搜索参数单元集P等一系列参数空间累积过程来确定d1与d3是否同圆，省去大量的无效累积过程[8-9]。

图3 利用梯度方向信息的改进RHT示例

另一方面，由于改进的RHT引入梯度方向信息进行判断，随机抽样点的个数从3个减少为2个，不考虑非圆点n，可知两点落在同一个待识别圆上的概率为：

可以看出，相对于3点，保证2点共圆的概率大幅增加，对于存在大量待识别圆的检测环境，这种概率的增加更为明显，无效累积现象明显减少。

1.3.2 利用先验知识减少无效累积

空间线性度检测层面上除了待识别的圆孔，还包含1个体模轮廓圆以及4个支撑柱圆，为剔除无效信息对圆检测过程的干扰，利用基于种子区域生长法对体模中心感兴趣区域进行提取，实现感兴趣区域的分割。分割过程如下：对图像进行二值化预处理后，在感兴趣区域设置初始种子像素点（图4）；通过检查初始种子点相邻像素的灰度特征，考虑是否将该相邻像素添加至生长区域，若种子与相邻像素的灰度差大于设定阈值，区域停止生长，得到体模中心感兴趣区域（图5）[10]。

由于感兴趣区域中还包含条状物模块，圆检测过程中也将会对其边缘点集进行提取。为避免累积无效圆参数，下面将利用半径的判断条件决定是否对圆参数p进行保留。

d1(x1,y1)为圆边缘点，d2(x2,y2)为矩形边缘点，两点梯度直线为l1、l2，两点连线中垂线为l3，此时l1、l2、l3交于点c(xc, yc)，若只考虑梯度信息，则认为两点可能同圆，并将候选圆参数p加入p，造成误判（图5）。本文将引入半径预设阈值o1，若边缘点d1(x1,y1)、d2(x2,y2)与交点c(xc, yc)的距离的距离大于预设阈值o1，则认为本次所选两点不共圆，无需对该圆参数进行累积，以提高检测效率（图6）。

图4 设置初始种子像素点

图5 种子生长结束所分割的感兴趣区域

图6 梯度信息的误判示例

综合上述改进措施，本文提出改进的RHT（RHT+）圆检测算法具体步骤，见图7。

图7 改进的RHT（RHT+）圆检测流程

1.4 空间线性度自动检测程序的实现

基于随机Hough变换圆检测算法的MRI QC空间线性度自动检测程序设计流程如下：

（1）对体模第四层横断位图像进行预处理。利用Otsu算法对图像进行二值化处理；针对图像上出现的噪点，可采用数学形态学方法（即选取合适尺度的结构元素，对图像进行膨胀、腐蚀操作）对其进行消除[11]。

（2）提取图像中模块的边缘信息。边缘检测算法利用边缘增强算子对边缘信息进行强化，并对像素边缘强度进行计算，最后通过与阈值的比较，实现对边缘点集的获取。本文采用边缘检测Canny算子进行边缘检测：应用一阶差分卷积模板实现目标边缘的增强；应用非极大值抑制实现局部梯度方向最大值的保留；应用双线性阈值实现边缘的检测和连接[12]。

（3）分别利用RHT与RHT+圆检测算法对目标圆孔进行检测，并对识别到的圆孔进行标记。

（4）已知目标圆孔位置，根据图片尺寸度量由式(1)分别计算频率与相位编码方向上的空间线性度。

2 结果

选择某次质量控制得到的Magphan SMR 170体模第二层图像作为检测对象，图像大小为256×256。算法编程环境为：MATLAB 2012b；运行硬件环境为：CPU 2.3 GHz、内存2 GB。

RHT与RHT+圆检测算法对目标圆孔的识别结果（识别到的圆孔用红圈标出），见图8～9；RHT+与RHT圆检测效果进行比较，见表1；人工（由熟练的质量控制工程师执行此次空间线性度检测程序）以及基于RHT+圆检测算法的线性度测量结果，见表2。

图8 RHT算法圆检测效果

图9 RHT+算法圆检测效果

表1 RHT+与RHT圆检测效果比较

综合上述实验结果：RHT与RHT+圆检测算法都实现了对13个目标圆孔的识别，相较于RHT算法，RHT+由于加入剔除无效累积的改进手段，检测速度获得近5倍的提升。表2中，基于RHT+算法的自动检测程序给出了空间线性度的测量结果，但该结果与人工检测结果存在一定出入，为比较二者差异，本文采用统计学方法行差异显著性检验。这里提出H0假设：自动检测结果与人工检测没有差异，显著性水平α设置为0.05级的情况下，P＞0.05时接受原假设，P＜0.05时拒绝原假设；检验结果得到P=0.10009＞0.05，接受H0假设，即认为两种方法检测结果差异无统计学意义[13]。

3 讨论

近几年，我国的MRI质量控制相关工作与研究呈现蓬勃发展之势，当中对于磁共振质量控制参数自动检测的讨论也愈发深入。作者官能成等[14]2015年在杂志发表质量控制自动检测相关研究成果，可以说是开创该研究领域之先河；直至2017年，张春青等[15]、李成等[16]重审研究背景，各自阐述了对质量控制自动检测系统构建方法的理解。以上等人的研究方向主要着眼于系统框架设计，未对参数检测所涉及的图像处理算法做具体讨论。而针对空间线性度，李馥然[17]利用边缘检测以及连通域标记等图像处理方法来实现该参数的自动检测，孟奥等[18]则利用人工手动选取目标区域法对线性度检测模块坐标进行定位。相较于以上研究方法，本文针对空间线性度提出了具体的检测算法，着力于检测过程人工的完全解放，在保证参数检测准确度的同时，强调检测效率，这对于后期构建全参数、实时的磁共振质量控制自动检测系统具有重要意义。

前文结果中所行差异显著性检验表明，两种方法检测结果差异无统计学意义，可以认为本文所设计的自动检测程序能够取代人工检测方式进行空间线性度的检测，这对于MRI质量控制程序的实施具备以下意义：① 摒弃了人工过程干预，无须受限于质量控制检测人员的检测行为习惯、操作熟练程度等主观因素影响，在追求随机误差最小化的过程中具备天然优势；② 算法运行的稳定性与连续性满足了质量控制程序对于检测连续性的要求，避免因为检测人员发生更替等原因而引起连续性的破坏，保证参数检测结果客观反映MRI系统相关性能状态的真实波动，这对于发现、隔离、解决成像系统细微问题具有重大意义；③ 规避质控程序对后处理工作站长时间占用的现象，解决系统校验行为与临床使用需求的冲突。

表2 人工与基于RHT+圆检测算法的线性度测量结果（%）

对于算法结构，本文所提出的空间线性度自动检测方法还应用了包括图像预处理、图像增强等算法，由于篇幅原因，本文对上述图像处理技术未做详细讨论，但不能否认其对自动检测目标实现的重要作用。此外，本文研究虽取得一定的成果，但若在更为复杂、严苛际成像环境中，如何在保证参数检测准确性的同时，增强算法的稳定性、鲁棒性、快速性，这是本文后续的研究目标，也是其他质量控制相关参数要实现自动检测的设计理念。