2. 广东省医学图像处理重点实验室,广东 广州 510515
2. Guangdong Provincial Key Laboratory of Medical Image Processing, Southern Medical University, Guangzhou 510515, China
磁共振成像(MR)能够较好的显示前列腺内在组织结构,对于前列腺疾病的分析与诊断具有重要的临床意义[1] 。临床上,对MR图像中前列腺的准确分割是对前列腺疾病进行诊断和病理阶段分析的重要前提,在前列腺切除术、放射治疗中也起着关键性的指导作用。由于前列腺形状不规则,边缘模糊,难以与周围其余组织区分开来,目前基于计算机的全自动分割算法误分割的概率大,离精确分割仍存在一定的差距。
目前前列腺MRI图像的分割算法主要有基于活动形状模型、活动表观模型、像素分类、多图谱等[2] 。基于形变模型[3-4] 的分割算法对初始形状的位置敏感,易陷入局部极值,且该算法的鲁棒性差、抗干扰性差。像素的分割算法很大程度上依赖于分类器和提取的特征的性能,分割效果差。基于多图谱[5] 的分割方法利用手动分割精度高的优势,将分割问题转化为配准问题。通过配准技术,它能有效的将手动分割好的医学图谱的形状先验知识融入分割过程,很好的重现分割结果。但当目标图像与图谱图像差异较大时,往往会因为配准误差过大而引入错误的形状先验知识信息,导致完全错误分割的情况[5-9] 。
为进一步校正配准引起的误差,实现准确、可靠的前列腺MR图像分割,受文献[10-11] 启发,本文将临床上医生的实时医学先验知识信息与多图谱分割算法相结合,提出了一种椭球先验约束下的前列腺MR图像多图谱分割算法。该椭球先验包含了医生或专家针对前列腺目标图像的形状和位置知识信息,可以有效的对通过配准引入的图谱形状先验知识进行校正和补偿。本文算法中,该形状先验的贡献主要表现在以下两个方面:(1)基于该椭球先验,可以只针对前列腺区域的相似度对图谱进行选择;(2)将该椭球先验引入图谱融合过程,对前列腺形状进行约束,使分割结果不再仅依赖于通过配准技术引入的前列腺图谱形状先验,有效确保了分割精度,同时也使得分割结果更加平滑。
1 椭球先验约束的前列腺MR图像分割算法 1.1 方法概述本文在以前提出的基于距离场融合的多图谱分割算法[12] 的基础上,引入了一种椭球形状先验对整个分割过程加以约束和校正。该算法主要可以分为以下六个步骤:预处理、图像配准、椭球形状先验生成、椭球先验约束的图谱选择、椭球先验约束的距离场融合、阈值分割等。图 1为本文方法的流程示意图。
为了避免偏移场、灰度对配准和分割的影响,实验中我们利用N3 算法[13] 对MR图像中的偏移场进行校正;并把图像灰度归一化到[0 100]。
1.3 图谱配准在基于多图谱分割的算法中,配准算法有:基于B样条的等网格配准方法,基于互信息的配准方法等。其中,基于互信息的配准方法有不依赖于图像本身灰度, 可实现自动校准、鲁棒性好等优点。本文中利用Klein 等人基于互信息的相似性测度开发出的配准软件ELASTIX[14] ,将前列腺图谱图像Ii与目标图像进行配准,并利用配准得到的空间变换T对对应的图谱前列腺标号图像Li进行映射,得到变形后的图谱标号图像Li '。
1.4 椭球形状先验的生成基于配准的多图谱分割算法的主要优势就是可以通过配准技术将专家或医生手动勾画好的医学图谱的形状先验知识和位置信息映射到目标图像上,从而实现目标图像的全自动分割。但当配准效果不好时,引入的这种形状先验信息往往是不准确的,导致分割精度大大降低。为了对不正确的前列腺形状先验知识信息进行补偿和校正,本文借助临床医生诊断时的医学先验知识信息引入了一种椭球形状先验。具体实现方法如下:
(1)对于每一个待分割目标图像,都由医生手工勾画出前列腺中间层上的四个边界点以及前列腺底部和尖端的中心点,其坐标可以依次表示(x1, y1, z1), ((x2, y2, z2)......(x6, y6, z6)。分别计算这些人工标记点的坐标在x,y,z方向的最小值和最大值,根据各个方向上的最小值和最大值的中点和距离长度则可以分别粗略估算出椭球的中心像素点坐标 (x0, y0, z0)以及三个方向上的轴半径 xr,yr,zr,利用公式(1)则可以初步确定前列腺的区域D。
$D = \left\{ {\left( {x,y,z} \right):\frac{{{{\left( {x - {x_0}} \right)}^2}}}{{x_r^2}} + \frac{{{{\left( {y - {y_0}} \right)}^2}}}{{y_r^2}} + \frac{{{{\left( {z - {z_0}} \right)}^2}}}{{z_r^2}} \le 1} \right\}$ | (1) |
其中,(x, y, z)为前列腺图像上某一像素点的坐标。
(2)对前列腺中间层、底部以及尖端的中心点坐标进行线性插值得到前列腺区域D 每层的中心点坐标(xi, yi, zi),并将每层的前列腺区域都分别对齐到该层的中心点 (xi, yi, zi)。
图 2为采用椭球拟合出的前列腺与专家手动勾画的前列腺的形状对比图。其中红色为专家手动分割结果,绿色为用椭球拟合出的结果。从图中可以看出,拟合出的椭球覆盖了大部分的前列腺区域。将该形状先验作为初始的分割结果,对前列腺区域加以限制和约束,可以很大程度上避免由配准误差所造成的完全错误分割,确保分割精度。
图谱选择就是从一系列图谱图像中,选出与目标图像最相像的图谱图像。2009年,Aljabar等[15] 通过实验证明图谱选择在基于多图谱的分割中起着至关重要的作用。尤其是在配准效果不太理想的时候,通过图谱选择可以剔除配准效果较差的图谱,很大程度上避免了配准误差对分割精度的影响。传统的多图谱分割算法中,图谱的选择多是基于全局图像的,但在前列腺MR图像中,前列腺区域只占图像的一小部分,若基于全局进行图谱选择,前列腺周围组织与器官会对图谱选择造成严重干扰,导致图谱选择不正确。如图 3A所示,前列腺区域(第2行)差异很大,但却由于周围组织相似,而被误选入图谱中。图 3B中,前列腺区域相似,但却由于周围组织相差较大,而漏选。
因此,本文算法中,不再是利用全局图像,而是针对前列腺区域进行图谱选择。配准后的图谱图像的前列腺区域可以由变形后的标号图像Li '得到,而目标图像的前列腺区域可以通过1.4 中生成的椭球先验D 得到。2000年,Artaechevarria等[16] 通过实验表明基于互信息的方法最能表达图谱选择中的近邻关系。因此本文算法采用Klein等[6] 提出的基于归一化互信息的相似性度量对前列腺区域进行相似性计算,从中选出互信息值不小于最大互信息值85%的图谱。
1.6 椭球先验约束的距离场融合算法标号融合算法就是对所选择的图谱进行线性加权组合的过程,标号融合算法的选取直接影响着最后的分割精度。
1.6.1 距离场融合算法在以前的工作中,我们提出了一种基于距离场的非局部块加权融合算法[13] (DistanceField Fusion, DFF)。基于距离场的非局部块加权标号融合算法(DFF)不再是对标号图像进行融合,而是对每个标号图像进行距离场变换后对相应的距离场图像进行融合。与目前国际上比较流行的标号融合算法(Major Voting [17], Weight Voting [18], SIMPLE [19], STAPLE[20], and Spatial STAPLE[21] 等)相比,它不再是融合单个像素点的标号,而是假设MR图像块构成的流形和距离场图像块构成的流形是一个局部微分同胚映射,由这个映射预测出每个MR图像块所对应的距离场图像块,再通过加权平均和阈值处理,最终得到MR图像的每个点所对应的标签。该算法的基本思想是:
(1)假设图谱图像块与其对应的距离场图像块分别位于两个非线性流形上,单个图像块可由其所在流形的近邻线性表示。依据假设,给定一个训练集
$\begin{array}{l} x_{test}^{DF}\left( p \right) = \sum\limits_{i = 1}^k {x_{ip}^{MR}} {w_{ip}} + \varepsilon \\ s.t\left\| \varepsilon \right\| \le \tau \\ \forall x_{ip}^{MR} \notin {N_k}\left( {x_{test}^{DF}\left( p \right)} \right),{\omega _{ip}} = 0 \end{array}$ | (2) |
其中,
(2)假设在局部区域约束条件下,这两个流形之间的映射近似为微分同胚映射。依据假设,在已知k近邻图谱图像块
$\begin{array}{l} x_{test}^{DF}\left( p \right) = f\left( {x_{test}^{MR}\left( p \right)} \right) \approx f\left( {\sum\limits_{i = 1}^k {x_{ip}^{MR}} {w_{ip}}} \right) = \\ \sum\limits_{i = 1}^k {f\left( {x_{ip}^{MR}} \right)} {w_{ip}} = \sum\limits_{i = 1}^k {x_{ip}^{DF}} {w_{ip}} \end{array}$ | (3) |
逐块预测待分割目标图像的距离场图像块,对于MR图像上的点p,其邻域上的点都会预测到一个距离场图像块,这些距离场图像块在p点位置上会有重叠,通过加权平均这些在p点位置重叠的距离场的值,即可得到p点的距离场,对每个点都进行以上相同的操作,即可得到整个目标图像的距离场图像F,再经过阈值分割,即可得到目标图像的标号图像。
1.6.2 椭球先验约束下的距离场融合算法基于距离场融合的多图谱分割算法引入了目标图像块的邻域信息和图谱标号距离场图像提供的位置信息,这种方法的优势就在于它降低了对图像配准的要求,即使配准存在一定的误差,也能通过搜索窗内的邻域信息以及图谱标号距离场图像提供的位置信息进行弥补。但由于个体之间前列腺大小、形状的差异性以及图谱的有限性,往往会出现配准误差过大的情况。当配准误差较大时,其分割结果与专家手动分割的结果仍然存在较大差异。图 7中可以看到配准效果不好时,基于距离场融合算法的部分较差的分割结果。
针对这个问题,我们将椭球先验作为初始分割结果引入图谱距离场融合过程,对图谱融合过程进行约束和补偿。加入了补偿项的距离场块预测方程式为:
$x_{test}^{DF}\left( p \right) = \sum\limits_{i = 1}^k {x_{ip}^{DF}} {w_{ip}} + \lambda \cdot {w_r}\left( p \right) \cdot d\left( p \right)$ | (4) |
其中d(p)为椭球先验对应的距离场图像中像素点p处对应的值,其定义如下:
$d\left( p \right) = \left\{ {\begin{array}{*{20}{l}} { - dis\left( {p,B} \right),poutsideC}\\ {dis\left( {p,B} \right),o{\rm{th}}erwise} \end{array}} \right.$ | (5) |
其中,C是椭球先验的曲面边界,B是C上距离点p最近的点,-dist(p, B)表示点p与点B的欧式距离。
λ为补偿项的全局权重。其表达式如下:
$\lambda = \exp \left( {1 - \frac{{2 \times \left| {CM \cap D} \right|}}{{\left| {CM \cap D} \right|}}} \right)$ | (6) |
其中D为1.4中椭球拟合的前列腺区域,
wr(p)为补偿项的局部权重,其定义如下:
${w_r}\left( p \right) = \left\{ {\begin{array}{*{20}{l}} {\frac{{d\left( p \right)}}{{{d_{\max }}}},0<d\left( p \right) \le {d_{\max }}}\\ { + \infty ,d\left( p \right) > {d_{\max }}}\\ {\frac{{d\left( p \right)}}{{{d_{\min }}}},{d_{\min }} \le d\left( p \right) \le 0}\\ { + \infty ,d\left( p \right)<{d_{\min }}} \end{array}} \right.$ | (7) |
其中dmax = pf ×Dmax;dmin = pb ×Dmin;Dmax, Dmin分别为椭球先验对应的距离场图像中的最大值、最小值;pf,pb为0-1之间的数。通过局部权重wr(p),则可以将待预测区域约束在椭球先验的一个窄带区域范围内,如图 4 所示。其中红色为专家手动勾画的结果,蓝色实线包围的区域为待预测区域。从图中可以看出,在局部权重的作用下,有效限制了前列腺区域范围,确保了分割精度。并且越靠近dmin所在边界的像素点补偿项的权重越大,即为非前列腺的概率越大,同理,越靠近dmax所在边界的像素点为前列腺的概率越大,从而在一定程度上保持了分割结果的连续性,也使得分割结果能在图谱融合项的作用下往真实边界靠拢。同时,也大大减少了要预测的像素点数,缩短了算法时间。
对公式(4)得到的距离场图像进行阈值处理,即可得到最终的分割结果L。
$L\left( p \right) = \left\{ {\begin{array}{*{20}{l}} {1,} & {F\left( p \right) > 0}\\ {0,} & {F\left( p \right) \le 0} \end{array}} \right.$ | (8) |
本文实验数据采用图像数据库MICCAI 2012公开的前列腺数据,该数据集包含50组T2加权的前列腺磁共振图像,每组数据包含图谱灰度图像及对应的专家手动勾画的标号图像。所有的实验均在MATLAB2014 中完成。实验中,采用leave-one-out 的方式对50 组前列腺数据进行分割。本文算法中的参数可以利用交叉验证法[21] 得到。本文中参数分别设置为:块大小5×5×3,搜索窗大小9×9×5,近邻k,pf ,pb分别取30,0.3,0.005。
图 5 是部分实验结果,从左往右依次为前列腺尖端、中间以及底部层的分割结果的截面图。其中红色为专家手动分割结果,蓝色为椭球拟合的初始分割结果,绿色为本文算法的分割结果。从图中可以看出,本文引入的椭球先验有效的确定了前列腺的位置,在椭球先验的基础上,本文算法可以有效的将前列腺区域分割出来。
为了验证本文算法图谱选择的有效性,图 6给出了本文的图谱选择方法与基于全局图像的图谱选择方法的图谱选择结果。从图中可以看出,与传统的基于全局图像选择的图谱图像相比,本文方法选择的图谱的前列腺区域与目标图像的前列腺区域更相似,表明了本文基于椭球约束区域进行图谱选择的有效性;同时,为了验证本文算法图谱融合过程中引入椭球先验约束的有效性,我们与之前提出的算法DFF进行比较。图 7为本文算法与DFF算法的最终分割结果。其中红色为专家手动分割结果,黄色为算法DFF的分割结果,绿色为本文算法的分割结果。从图中可以看出,本文算法在椭球形状的约束下分割结果与前列腺的真实边界更加接近,特别是在目标图像边界不清楚的情况下,本文算法通过椭球先验引入的形状和位置信息,成功的将演化曲线吸引到了真实的前列腺边界。同时,本文算法还有效改善了DFF算法分割结果在边界会出现一些毛刺或凹陷现象的情况,可以得到更加光滑、连续的边界。
为了对分割算法的准确程度有一个客观的评价,本文将分割结果与专家勾画结果进行比较。采用Dice相似性系数(Dice Similarity Coefficient,DSC)作为衡量标准,其定义为:
$DSC\left( {A,B} \right) = \frac{{2 \times \left( {A \cap B} \right)}}{{A \cup B}}$ | (9) |
其中A和B分别是算法分割结果和专家手动勾画结果,A ∩ B表示算法分割结果和专家手动勾画结果重叠部分的体积,A ∩ B表示算法分割结果和专家手动勾画结果体积之和。图 8给出了本文算法分割结果和DFF算法分割结果分别与专家手动勾画结果的Dice 相似性系数。从图中可以看出,用本文算法分割的50例数据的Dice 相似性系数均在80%以上,且平均精度可达到88.12%;基于距离场融合的算法在不同个体上的分割精度上下浮动较大,平均精度只有83.01%;而其他基于多图谱的前列腺MR图像分割算法的文献[5-7] 中,前列腺平均分割精度不高,介于0.75~0.8之间;文献[8-9] 提出的基于多图谱的前列腺分割算法中,平均分割精度虽有所提高(84%左右),但是仍然存在部分错误分割(DSC<0.2)的情况。综上可以看出,本文算法有效的避免了由配准误差过大引起的错误分割的情况,具有较高的分割精度以及较好的一致性。表 1列出本文方法和目前已公开的效果比较好的其他前列腺分割方法在本文数据集上的分割性能,本文方法对应的DSC指标优于其他方法。
最后,为了评估人工选取的标记点对算法的影响,实验中我们分别在前列腺底部,中间以及尖端层的上下相邻层上,随机选取1 层勾画人工标记点,进行实验。图 9为10次实验的结果。从图中可以看出,本文算法对人工标记点的选取并不敏感,具有较好的鲁棒性。
针对前列腺多图谱分割算法分割精度易受配准效果的影响,鲁棒性差的问题,本文将多图谱分割与椭球形状先验知识相结合,提出了一种椭球先验约束下的前列腺MR图像分割算法。与其他基于多图谱的前列腺MR图像分割算法相比,本文算法有效的避免了由配准误差引起的错误分割的情况,具有较好的一致性,且明显提高了分割精度。由此表明本文提出的算法是一种有效的前列腺MR图像分割算法。后续工作中,可将本文算法与其他图像分割算法,如图割等相结合,进一步提高分割精度。
[1] | 曹一挥. 基于多图谱医学图像分割的技术研究[D].西安: 中国科学院大学, 2013. |
[2] | Litjens G, Toth R, van de Ven W, et al. Evaluation of prostate segmentation algorithms for MRI: the PROMISE12 challenge[J]. Med Image Anal, 2014, 18 (2): 359-73. DOI: 10.1016/j.media.2013.12.002. |
[3] | Kirschner M, Jung F, Wesarg S. Automatic prostate segmentation in Mr images with a probabilistic active shape model: the PROMISE12 challenge[EB/OL] (2012-06-29) [2015-05-14]. https://grand-challenge.org/site/promise12/serve/public_html/pdfs/grislies.pdf. |
[4] | Maan B. Prostate MR image segmentation using 3D active appearance models: the PROMISE12 challenge [EB/OL]. (2012-06-29) [2015-05-14]. https://grand-challenge.org/site/promise12/serve/public_html/pdfs/Utwente_Signals.pdf. |
[5] | Klein S, van der Heide UA, Lips IM, et al. Automatic segmentation of the prostate in 3D Mr images by Atlas matching using localized mutual information[J]. Med Phys, 2008, 35 (4): 1407-17. DOI: 10.1118/1.2842076. |
[6] | Martin S, Daanen V, Troccaz J. Atlas-based prostate segmentation using an hybrid registration[J]. Int J Comput Assist Radiol Surg, 2008, 3 (6): 485-92. DOI: 10.1007/s11548-008-0247-0. |
[7] | Dowling J, Fripp J, Greer P, et al. Automatic atlas-based segmentation of the prostate:a miccai 2009 prostate segmentation challenge entry[EB/OL] [2017-03-06].http://www.wiki.na-mic.org/Wiki/images/f/f1/Dowling_2009_MICCAIProstate_v2.pdf. |
[8] | Ou, Y, Doshi, et al. Multi-Atlas segmentation of the prostate:a zooming process with robust registration and Atlas selection:the PROMISE12 challenge[EB/OL].(2012-04-07) [2015-05-14]. https://grand-challenge.org/site/promise12/serve/public_html/pdfs/SBIA.pdf. |
[9] | Gao, Q, Rueckert, et al. An automatic multi-atlas based prostate segmentation using local appearance-specific atlases and patchbased voxel weighting: the PROMISE12challenge [EB/OL].(2012-02-07) [2015-05-14]. https://grand-challenge.org/site/promise12/serve/public_html/pdfs/ICProstateSeg.pdf. |
[10] | Slabaugh G, Unal G. Graph cuts segmentation using an elliptical shape prior[C] //2005 international conference on image processing (ICIP), VOLS 1-5. 2, 2005: 1973-6. |
[11] | Khalvati F, Salmanpour A, Rahnamayan S, et al. Sequential Registration-Based segmentation of the prostate gland in Mr image volumes[J]. J Digit Imaging, 2016, 29 (2): 254-63. DOI: 10.1007/s10278-015-9844-y. |
[12] | Pang S, Lu Z, Yang W, et al. Hippocampus segmentation through distance field fusion[M]. 2015: 104-11. |
[13] | Nyúl LG, Udupa JK. On standardizing the Mr image intensity scale[J]. Magn Reson Med, 1999, 42 (6): 1072-81. DOI: 10.1002/(ISSN)1522-2594. |
[14] | Klein S, Staring M, Murphy K, et al. Elastix: a toolbox for intensitybased medical image registration[J]. IEEE Trans Med Imaging, 2010, 29 (1): 196-205. DOI: 10.1109/TMI.2009.2035616. |
[15] | Aljabar P, Heckemann RA, Hammers A, et al. Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[J]. Neuroimage, 2009, 46 (3): 726-38. DOI: 10.1016/j.neuroimage.2009.02.018. |
[16] | Artaechevarria X, Munoz-Barrutia A, Ortiz-de-Solorzano C. Combination strategies in multi-atlas image segmentation: application to brain Mr data[J]. IEEE Trans Med Imaging, 2009, 28 (8): 1266-77. DOI: 10.1109/TMI.2009.2014372. |
[17] | Cabezas M, Oliver A, Lladó X, et al. A review of atlas-based segmentation for magnetic resonance brain images[J]. Comput Methods Programs Biomed, 2011, 104 (3): e158-77. DOI: 10.1016/j.cmpb.2011.07.015. |
[18] | Wu G, Wang Q, Zhang D, et al. A generative probability model of joint label fusion for multi-atlas based brain segmentation[J]. Med Image Anal, 2014, 18 (6): 881-90. DOI: 10.1016/j.media.2013.10.013. |
[19] | Langerak TR, van der Heide UA, Kotte AN, et al. Label fusion in atlas-based segmentation using a selective and iterative method for performance level estimation (SIMPLE)[J]. IEEE Trans Med Imaging, 2010, 29 (12): 2000-8. DOI: 10.1109/TMI.2010.2057442. |
[20] | Warfield SK, Zou KH, Wells WM. Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation[J]. IEEE Trans Med Imaging, 2004, 23 (7): 903-21. DOI: 10.1109/TMI.2004.828354. |
[21] | Asman AJ, Landman BA. Formulating spatially varying performance in the statistical fusion framework[J]. IEEE Trans Med Imaging, 2012, 31 (6): 1326-36. DOI: 10.1109/TMI.2012.2190992. |