本发明涉及医学图像处理技术领域，具体涉及一种基于凸约束种子区域生长的脊髓分割方法。
背景技术：
随着医学影像在临床医学诊断中的广泛应用，图像分割技术作为医学图像处理中的关键一步，其主要目的是对感兴趣目标如脊髓、脑组织、血管、细胞等进行准确的分割，在医学中发挥着越来越大的作用。医学图像分割是提取医学影像中特殊组织定量信息所不可缺少的手段，其被广泛地应用于各个方面，如组织容积的定量分析、解剖结构的研究、病变组织的定位、治疗规划、功能成像数据的局部体效应校正、计算机辅助诊断、计算机引导手术等等。
在医学图像分割领域，比较经典的分割算法有阈值分割方法、水平集方法、活动轮廓模型、基于图论的分割方法、基于智能群优化的分割方法、基于深度学习的分割方法、种子区域生长等。在ct图像对脊髓的分割中，请参考图2所示，可以发现待分割脊髓区域在亮度上与其他区域非常相近，因此阈值分割方法和水平集方法并不适合处理此类图像。活动轮廓模型是利用梯度信息进行曲线演化来完成区域分割，但从图2中可以发现，脊髓区域边界并不十分明显，因此采用活动轮廓模型很容易导致边缘泄露。而基于图论的分割方法是以图像中所有像素为节点建立图模型进行分割，实时性较差，并且在图2中很难建立较好的前景和背景模型。而基于智能群优化的分割方法需要大量的迭代次数寻求最优解，同样实时性较差。基于深度学习的分割方法是目前的一个研究热点，但是该方法需要大量的样本图像进行训练，并且每幅图像均需要人工进行标记，因而对训练样本要求较高且成本较大。观察图2可以发现，虽然脊髓区域的亮度与其他区域亮度具有很高的相似性，但与其周围却有一定的差异，如果能仅对这一局部区域进行处理便可有效完成脊髓区域的提取。因此，针对脊髓分割问题种子区域生长方法将是一个不错的选择方案。
本发明的发明人经过研究发现，传统的种子区域生长方法是仅利用图像灰度信息，每次优先生长与种子区域相邻且灰度值最相近的候选点，而没有充分利用其它信息，其生长过程无规律可寻且不集中，有时会沿着某一条路径一直延伸下去，从而导致在脊髓区域内部还没有生长完全时就发生了过生长(即生长泄露)现象。因此，针对传统种子区域生长方法在生长过程中如何保持集中生长并避免过早出现过生长现象，以及阻止过生长的发生，成为目前种子区域生长方法急需解决的问题。
技术实现要素：
针对现有种子区域生长方法在生长过程中如何保持集中生长并避免过早出现过生长现象，以及阻止过生长发生的技术问题，本发明提供一种基于凸约束种子区域生长的脊髓分割方法，该方法通过使种子区域在生长过程中尽可能保持凸性质，从而有效地完成了ct图像中脊髓区域的准确分割。
为了解决上述技术问题，本发明采用了如下的技术方案：
一种基于凸约束种子区域生长的脊髓分割方法，包括以下步骤：
s1、对输入ct图像进行对比度增强，以提高图像中脊髓区域与其周围区域在亮度上的差异，并通过计算图像空间位置质心来获取初始种子点；
s2、在种子区域生长过程中，选择具有最大曲率且与种子区域在灰度信息上具有最大相似度的候选种子点作为下一次的生长点，即优先生长那些处于凹形区域内的候选种子点；
s3、利用种子区域的区域曲率来表示种子区域的形状信息，通过检测种子区域形状的变化来判断是否发生了生长泄露，若发生了生长泄露则停止生长并提取脊髓区域，若未发生生长泄露则返回步骤s2继续生长。
进一步，所述步骤s1中，对输入ct图像进行对比度增强包括：
将ct图像数据通过下式归一化处理到[0，1]：
其中，src为输入图像，nor为归一化后图像，i为任意一像素点，min(src)为输入图像src的最小值，max(src)为输入图像src的最大值；
对归一化后图像采用下式进行对比度调整：
其中，i为对比度调整后的输出图像，τ为对比度调整阈值，且τ∈[0,1]。
进一步，所述步骤s1中，图像空间位置质心的计算如下：
其中，(x(i)，y(i))为像素点i的位置坐标，(cx，cy)为图像i的图像空间位置质心。
进一步，所述步骤s2中，采用曲率cur(c)来度量某一候选种子点c的凸凹程度，定义如下：
其中，
进一步，所述步骤s2中，下一次生长的候选种子点通过下式获得：
其中，
进一步，所述步骤s3中，优选利用候选种子点
其中，|ri|表示种子点个数，
进一步，所述步骤s3中，通过检测种子区域形状的变化来判断是否发生了生长泄露包括：
s31、判断种子区域内的种子数量是否大于预设阈值nδ，若为是则执行步骤s32，否则返回步骤s2继续生长；
s32、计算第i次生长结果得到的种子区域ri的区域曲率con(ri)，若con(ri)小于生长过程中的最小区域曲率mincon，则更新mincon＝con(ri)；再判断是否满足con(ri)-mincon≥δ，其中δ为预设的最大容忍误差，若当前区域曲率与最小区域曲率的差超过最大容忍误差则停止生长，否则返回步骤s2继续生长。
与现有技术相比，本发明提供的基于凸约束种子区域生长的脊髓分割方法，通过使种子区域在生长过程中尽可能保持凸性质，从而能更为有效地完成ct图像中的脊髓分割，具体在种子区域生长过程中，通过加入凸约束来选择具有最大曲率且与种子区域在灰度信息上具有最大相似度的候选种子点作为下一次的生长点，有效实现了种子区域在生长过程中保持集中生长并避免了过早出现过生长现象，使得在生长过程中种子区域能够保持较好的凸形状；再者，随着种子生长的进行，当目标区域内部生长完全后一旦发生了生长泄露，就会造成种子区域的形状发生较大的变化，即种子区域由凸形状过渡成了凹形状，而本发明同样利用区域曲率来表示种子区域的形状信息，并通过检测种子区域形状的变化来判断是否出现了生长泄露，以此作为停止生长的依据，即若出现了生长泄露则停止生长，反之则继续生长，因而有效阻止了过生长发生。
附图说明
图1是本发明提供的基于凸约束种子区域生长的脊髓分割方法流程示意图。
图2是本发明提供的一待处理原始ct图像。
图3是对原始ct图像经过对比度调整后的图像。
图4是对原始ct图像经过对比度调整后的图像及其计算出的质心位置。
图5是ct图像f1_5的部分生长结果示意图。
图6是不同ct图像中脊髓区域的比较示意图。
图7是种子区域与候选种子区域示意图。
图8是不同候选种子点示意图。
图9是本方法与传统方法生长过程比较示意图。
图10是本方法与传统方法最终生长结果比较示意图。
图11是种子区域、边缘种子点以及候选种子点的示意图。
图12是不同ct图像种子区域的曲率变化结果示意图。
图13是ct图像f1_5的种子区域曲率局部变化示意图。
图14是ct图像f1_5的部分生长过程示意图。
具体实施方式
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体图示，进一步阐述本发明。
请参考图1至图14所示，本发明提供一种基于凸约束种子区域生长的脊髓分割方法，包括以下步骤：
s1、对输入ct图像进行对比度增强，以提高图像中脊髓区域与其周围区域在亮度上的差异，并通过计算图像空间位置质心来获取初始种子点。作为具体实施例，所述对输入ct图像进行对比度增强包括以下步骤：
将ct图像数据通过下式归一化处理到[0，1]区间，由此实现了在实验过程中方便统一处理：
其中，src为输入图像，nor为归一化后图像，i为任意一像素点，min(src)为输入图像src的最小值，max(src)为输入图像src的最大值；
图2为本发明提供的一幅待处理ct图像，从该图中可以发现，虽然脊髓部分与其周围区域在亮度上存在一定差异，但是在对比度上并不十分明显，因此本发明首先对待处理图像进行对比度增强，具体对归一化后图像采用下式进行对比度调整：
其中，i为对比度调整后的输出图像，τ为对比度调整阈值，且τ∈[0,1]；作为一种实施方式，所述对比度调整阈值τ优选设为0.5。如图3所示为对图2利用对比度调整阈值τ＝0.5进行对比度调整后的ct图像，从图3中可以发现，经过对比度调整后脊髓区域与其周围区域在亮度上差异更为明显，在此基础上即可进行有效的种子区域生长。
观察测试图像可以发现脊髓周围区域具有较大的亮度值，通过多次调整阈值τ其周围区域仍能得到较为完整的保留。如图4所示为当τ＝0.7时所得到的ct图像，观察此图可以发现，利用亮度加权得到的中心位置位于脊髓区域内部，根据这一发现，本发明通过计算图像空间位置质心来作为初始种子点，具体所述图像空间位置质心的计算如下：
其中，(x(i)，y(i))为像素点i的位置坐标，(cx，cy)为图像i的图像空间位置质心。
s2、利用ct图像中每个脊髓区域均保持一定的凸性这一特征，在种子区域生长过程中，通过加入凸约束来选择具有最大曲率且与种子区域在灰度信息上具有最大相似度的候选种子点作为下一次的生长点，即优先生长那些处于凹形区域内的候选种子点。
具体地，在采用传统种子区域生长算法进行区域生长时发现，其生长过程无规律可循，并且有时会沿某一条路径一直延伸下去，如图5所示为一测试ct图像f1_5的部分生长结果，其中(a)是生长200次的结果，(b)是生长230次的结果，(c)是生长250次的结果，(d)是是生长300次的结果。从该图中可以看出，其生长过程并不集中，另外在脊髓区域内部出现了还没有生长完全就导致了过生长(生长泄露)现象的发生。因此，本发明需要解决在生长过程中如何保持集中生长以及避免过早出现过生长两个问题。
不同ct图像中的脊髓区域如图6所示，其中(a)为ct图像f1_1及其脊髓区域，(b)为ct图像f1_2及其脊髓区域，(c)为ct图像f1_3及其脊髓区域，(d)为ct图像f1_4及其脊髓区域，(e)为ct图像f1_5及其脊髓区域，通过观察可以发现每个脊髓区域均保持一定的凸性。然而本发明的发明人经过研究发现，在传统种子区域生长算法的生长过程中均没有有效利用这种凸性。如果在生长过程中尽可能保持种子区域具有一定的凸性，那么便可达到集中生长以及避免过早出现生长泄露现象。
通过前述观察，本发明关键性地提出了在生长过程中保持种子区域具有一定的凸性，即优先生长那些处于凹形区域内的候选种子点。如图7所示为某次生长结果中的种子区域ri(黑色)与候选种子区域candidate(ri)(灰色)示意图，其中候选种子的获取采用了4邻域系统。
作为具体实施例，所述步骤s2中，采用曲率cur(c)来度量某一候选种子点c的凸凹程度，记为cur(c)，定义如下：
其中，
记像素点p的灰度值为i(p)，种子区域ri的灰度均值为mean(ri)，候选种子区域为
而本发明要求下一次生长的候选种子点还应具有较大的曲率，因此作为具体实施例，所述步骤s2中，下一次生长的候选种子点通过下式获得：
其中，
因此，与传统种子生长方法相比，本发明方法的生长性能得到了明显提升，如图10所示为最终生长结果的对比，其中(a)为传统方法生长结果，(b)为本发明方法生长结果。而本发明的发明人经过进一步研究发现，虽然前述方法能使生长过程比较集中，但仍无法避免过生长现象的发生。那么如何阻止过生长现象发生仍是本发明需要解决的另一个关键问题，对此本发明方法还包括步骤：s3、利用种子区域的区域曲率来表示种子区域的形状信息，通过检测种子区域形状的变化来判断是否发生了生长泄露，若发生了生长泄露则停止生长并提取脊髓区域，若未发生生长泄露则返回步骤s2继续生长。
从前述图10中的实验结果可以发现，种子生长过程比较集中始终保持了较好的凸性，并且当目标区域内部生长完全后才造成生长泄露，因此只需设置合适的停止生长时间来阻止进一步生长即可。但是，在前述步骤s2中仅考虑到候选种子点的曲率信息，而没有整体考虑到种子区域的形状信息。观察图10(b)可以发现，在生长泄露发生时会造成种子区域的形状发生较大的变化，因此只要充分利用形状变化这一有用信息即可有效阻止种子区域继续生长。
根据前述分析，本发明将使用种子区域ri的凹凸度con′(ri)(区域曲率)来表示种子区域的形状信息，即利用所有种子点记为
其中，|ri|表示种子点个数，cur(pij)表示种子点pij的曲率。记边缘种子点为
当某一边缘像素点
在生长过程中，随着种子生长的进行，其种子点数也逐渐增多，若按公式(9)计算每次生长结果的凹凸度，均需要先获取所有边缘种子点，其计算过程较为复杂。然而在生长过程中，由于每次生长得到的候选种子点均已保存在一个数组中，因此本发明将创造性地利用候选种子点来完成种子区域凹凸度即区域曲率的近似计算。记种子区域ri的候选种子点为
其中，|ri|表示种子点个数，
与con′(ri)相反，当con(ri)越小时说明种子区域ri越凸，反之种子区域ri越凹。随着种子生长的进行，区域会保持更好的凸性质，一旦生长泄露发生，种子区域的凸性将被破坏，以此便可检测出是否发生了生长泄露。如图12所示为五幅测试图像生长过程的种子区域曲率变化趋势，其中(a1)为f1_1生长结果，(b1)为f1_1种子区域曲率变化，(a2)为f1_2生长结果，(b2)为f1_2种子区域曲率变化，(a3)为f1_3生长结果，(b3)为f1_3种子区域曲率变化，(a4)为f1_4生长结果，(b4)为f1_4种子区域曲率变化，(a5)为f1_5生长结果，(b5)为f1_5种子区域曲率变化。从前四幅测试图像可以看出，在生长过程中种子区域曲率呈下降趋势减小，直到填充整个脊髓区域；而对于第五幅测试图像，在出现生长泄露前呈下降趋势，当发生泄露后，呈逐渐上升变化趋势。如图13所示为测试图像f1_5种子区域生长曲率的局部变化，其中(a)为局部变化曲线1，(b)为局部变化曲线2。从该图中可以看出：在生长230次与270次之间发生较为明显的跳变，并从生长270次之后开始呈上升趋势，过生长也随之开始。如图14所示，分别为测试图像f1_5生长250次、260次、270次、280次290次以及300次后的区域生长结果。
从图12前四幅测试图像的种子区域曲率变化可以发现，在初始生长阶段并不呈明显的下降趋势，稳定性较差，当达到一定生长次数后，曲率变化成明显下降趋势。另外，步骤s2中的策略避免了过早出现生长泄露，而直至脊髓区域生长完全后才会发生过生长。因此，在生长过程中只需当种子区域种子数量达到一定数量nδ后就开始检测是否出现过生长。在实际应用中nδ可根据经验或训练方式获取，作为一种实施方式nδ＝200。
作为具体实施例，所述步骤s3中，通过检测种子区域形状的变化来判断是否发生了生长泄露，即过生长检测过程包括：
s31、判断种子区域内的种子数量是否大于预设阈值nδ，若为是则执行步骤s32，否则返回步骤s2继续生长；
s32、计算第i次生长结果得到的种子区域ri的区域曲率con(ri)，若con(ri)小于生长过程中的最小区域曲率mincon，则更新mincon＝con(ri)；再判断是否满足con(ri)-mincon≥δ，其中δ为预设的最大容忍误差，若当前区域曲率与最小区域曲率的差超过最大容忍误差则停止生长，否则返回步骤s2继续生长；其中，最大容忍误差δ可根据经验或训练方式获取，本实验中取δ＝0.007。具体在对测试图像f1_5进行生长时，如图13中(a)所示，当生长到第229次时获得最小值，而当生长到第256次时由于超过了给定的最大容忍误差，因而生长停止。
与现有技术相比，本发明提供的基于凸约束种子区域生长的脊髓分割方法，通过使种子区域在生长过程中尽可能保持凸性质，从而能更为有效地完成ct图像中的脊髓分割，具体在种子区域生长过程中，通过加入凸约束来选择具有最大曲率且与种子区域在灰度信息上具有最大相似度的候选种子点作为下一次的生长点，有效实现了种子区域在生长过程中保持集中生长并避免了过早出现过生长现象，使得在生长过程中种子区域能够保持较好的凸形状；再者，随着种子生长的进行，当目标区域内部生长完全后一旦发生了生长泄露，就会造成种子区域的形状发生较大的变化，即种子区域由凸形状过渡成了凹形状，而本发明同样利用区域曲率来表示种子区域的形状信息，并通过检测种子区域形状的变化来判断是否出现了生长泄露，以此作为停止生长的依据，即若出现了生长泄露则停止生长，反之则继续生长，因而有效阻止了过生长发生。
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。