一种基于多角度序列化图像空间特征点集的脑血管图像分割方法
【技术领域】
[0001] 本发明涉及图像处理技术领域，尤其涉及一种基于多角度序列化图像空间特征点 集的脑血管图像分割方法。
【背景技术】
[0002] 通过DSA技术得到脑血管的图像，可以帮助医生对患者进行诊断，但是由于患者 的呼吸吞咽等动作，会存在明显的伪影和噪声，影响了医生的观测和诊断，因此，提取准确 的血管组织图像是当前基于DSA技术的血管诊疗中急需解决的问题。现有技术中包括手动 分割方法，配准方法，基于模型的方法。手动分割的方法需要人工进行交互，不同的人工操 作得到的分割结果差别很大，还需要耗费大量的时间，通过配准方法得到的血管图像中包 含较大的噪声，影响医生的诊断，采用基于模型的方法，需要一张模板图像，只有精确的模 板才能得到较好的血管分割结果。
【发明内容】
[0003] 根据现有技术存在的问题，本发明公开了一种基于多角度序列化图像空间特征点 集的脑血管图像分割方法，包括以下步骤：
[0004] Sl :对脑血管的活片图像和蒙片图像进行配准：在活片图像和蒙片图像上确定几 何特征点，找出匹配的特征点对，通过匹配的特征点对之间的关系对蒙片图像进行处理获 得序列化减影图像；
[0005] S2 :提取序列化减影图像的几何特征点集；
[0006] S3 :对序列化减影图像中血管边缘的特征点进行局部位置调整：通过灰度值梯度 来调整特征点的位置，将血管边缘的特征点移动到血管内；
[0007] S4 :采用空间旋转坐标系模型剔除减影图像中非血管位置的特征点：建立一个空 间旋转坐标系模型，利用相邻减影图像上特征点的位置信息，求出特征点在其它减影图像 中的理论位置，根据其他减影图像中特征点的位置，采取投票的方式剔除非血管位置的特 征点；
[0008] S5 :确定序列化减影图像的自适应分割阈值；
[0009] S6 :基于区域生长和自适应阈值的血管图像分割：将S4中获得的特征点作为种子 点，将S5中获得的自适应分割阈值作为生长法则的标准，采用区域生长算法分割序列化减 影图像获得纯净的脑血管图像。
[0010] Sl中具体采用如下方式：
[0011] Sll :利用SURF算法分别在相对应的活片图像和蒙片图像上确定几何特征点，计 算出特征点的特征向量；
[0012] S12 :利用FLANN匹配器进行特征向量匹配，得到若干匹配点对；
[0013] S13 :对匹配好的点对进行筛选消除错误匹配，计算出匹配点对之间的最小距离， 只保留距离小于匹配点对之间的最小距离三倍的点对；
[0014] S14 :采用RANSAC算法去除错误匹配点对，得到两幅图像之间的变换矩阵，通过变 换矩阵获得配准后的蒙片图像，将活片图像与配准后的蒙片图像进行直接减影，得到序列 化减影图像。
[0015] S2中具体采用如下方式：
[0016] S21 :利用SIFT算法在序列化减影图像上和蒙片图像上提取几何特征点集；
[0017] S22:在获得减影图像的特征点时对图像进行上下分片操作，取图像的三分之一作 为上片图像，其余作为下片图像，对上片图像进行增强对比度的操作获得较多的特征点；
[0018] S23:计算蒙片图像中的特征点与减影图像中的特征点之间的欧式距离，如果减影 图像中的特征点与蒙片图像中的特征点之间的欧式距离小于50,则将该特征点其剔除。
[0019] S3中具体采用如下方式：
[0020] S31 :判断减影图像中特征点集是否都被处理完毕，如果是则计算过程结束，如果 不是则从特征点集中取出一个未被处理的特征点；
[0021] S32 :判断上述未被处理的特征点是否为8领域内灰度值最小的点，如果是则继续 判断下一个未被处理的特征点；
[0022] S33 :如果S32中判断为否则计算出该特征点8个方向上的梯度，将梯度最大的点 作为新的特征点；
[0023] S34:重复上述步骤，直到减影图像中特征点集都处理完毕。
[0024] S4中具体采用如下方式：
[0025] S41 :从当前的减影图像中选取一特征点P1 (X1, yi)，在当前的减影图像中的相邻减 影图像中计算与P1满足如下关系的特征点P 2(x2, y2):
[0026] X2-X11 ^ 4
[0027] Y2-Y11 ^ 2
[0028] S42 :根据横坐标横坐标X2以及P JP P2之间的角度α = 2°计算特 征点P1到旋转中轴的距离r:
[0029]
[0030] S43:计算特征点P1在其他减影图像中投影的理论位置Pn:设当前处理的减影图像 A的序号为其他减影图像B的序号为、，通过序号的差值求出两幅减影图像之间的旋转 的角度α :
[0031] a = |2X (In-I1) | ；
[0032] S44 :根据上述计算出的r和旋转角度α求出特征点？1在当前减影图像B中的投 景多的位置Pn(XnJ n);
[0033] Xn=XiirXcosa
[0034] yn= y 1
[0035] S45 :从减影图像B中查找是否存在满足以下要求的特征点（Xni，yj，如果存在则给 特征点P1投一票，否则不投票；
[0036] xn-xnI ^ 4
[0037] yn-yn| ^ 2
[0038] S46 :设定阈值，剔除得票数低于该阈值的特征点后获得特征点集。
[0039] S5中具体采用如下方式：
[0040] S51 :对减影图像进行高斯滤波操作，去除噪声，采用如S22中所述的方式对图像 进行分片处理，计算上片图像的自适应分割阈值：
[0041] 以像素点P (x，y)为中心，取一个wupX Wup的窗口，其中Wup为奇数；
[0042] 计算窗口内的灰度期望值Elcical (X，y);
[0043]
[0044]
[0045]
[0046] 根据窗口内的灰度期望值Elocal (X，y)及窗口内灰度分布标准差〇 (X，y)，确定 阈值 Tup(x，y):
[0047] Tup (x, y) = Elocal (x, y) +kup X σ (χ, y)
[0048] 其中：kup取值范围为0. 1-0. 3, wup取值范围为6-11 ;
[0049] 计算下片图像的自适应分割阈值：
[0050] 计算下片图像的积分图，积分图中任一点I(x，y)的值等于从图像左上角到这个 点所构成的矩形区域内的所有的点的灰度值之和，I表示积分图像，G表示减影图像，其中 χ, O^j ^y〇
[0051] I (x, y) = sum(G(i, j))
[0052] 以像素点P (x，y)为中心，取一个wd_X wd_的窗口，其中w d_为奇数，
[0053] 通过积分图得到窗口内的像素的平均值m(x, y)，其中x-wd_/2彡i彡x+wd_/2, y-w d〇X j 彡 y+wdown/2〇
[0054]
[0055] m(x,y) = Sum(G(i, j))/wdown2
[0056] 根据窗口内的平均灰度值m(x, y)，确定阈值Td_ (x, y):
[0057] Tdown (x, y) = kdownXm(x, y)
[0058] 其中：kd_取值范围