CN105956587A - 一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法 - Google Patents

Abstract

一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，本发明涉及基于形状约束的膝关节磁共振图像序列半月板自动提取方法。本发明为了解决现有半月板分割方法计算效率与处理效果的不平衡、利用单一灰度信息使处理效果不理想、需要人工设置初始点、难以实现自动图像序列批处理、提取效果受成像条件影响等问题。具体过程为：一、ROI定位，选定初始片层组；二、进行优化，得到优化后片层集合；三、对优化后片层集合中的片层的ROI进行阈值自动选择与半月板目标提取；四、根据提取出的半月板目标在优化后片层集合中进行半月板原始信息提取，得到包含分割出的半月板原始信息的ROI的片层。本发明应用于半月板提取领域。

Description

一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法

技术领域

本发明涉及基于形状约束的膝关节磁共振图像序列半月板自动提取方法。

背景技术

现有的膝关节半月板提取方法有基于目标边界信息的的分割方法，包括边缘检测算子、区域生长、分水岭算法等。此类方法施行起来较为简单、运算量小，但利用的信息较单一、无法结合医学经验知识，目标分割效果差；有基于轮廓模型的方法，包括snake模型及其改进算法，如主动轮廓模型、测地线主动轮廓模型等。此类方法处理效果较好，但计算量大、模型演化速度慢、初始条件需要人工设置且对最终处理结果影响大，同时难以实现图像序列批处理；有基于灰度信息的提取方法，包括纹理目标分割、阈值处理等。此类方法应用了更多的图像特征及医学经验知识，可以较好的平衡计算效率与处理效果，但存在受成像条件影响产生误提取的现象，如ROI定位受成像位置影响、关节软骨被误提取为半月板目标等现象。

发明内容

本发明为了解决现有半月板分割方法中存在的以下问题：计算效率与处理效果的不平衡、未结合医学经验知识利用单一灰度信息使处理效果不理想、需要人工设置初始点、难以实现自动图像序列批处理、提取效果受成像条件影响等问题，从而提出一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、ROI(感兴趣区域)定位，选定初始片层组S_selected；

所述，ROI为感兴趣区域；

步骤二、对选定的初始片层组S_selected进行优化，得到优化后片层集合S_optimized；

步骤三、对优化后片层集合S_optimized中的片层的ROI进行阈值自动选择与半月板目标提取；

步骤四、根据提取出的半月板目标在优化后片层集合S_optimized中进行半月板原始信息提取，得到包含分割出的半月板原始信息的ROI的片层T_k。

发明效果

实验中采用的三组实验数据是由哈尔滨医科大学附属医院提供，采集于三位年龄分布在30-50之间、于2014年在哈尔滨医科大学附属医院进行膝关节磁共振检查的患者，图像为T1加权(TR＝500ms，TE＝17ms)，片层间隔为4mm，像素尺寸为512*512。为尊重用户隐私，成像序列分别称为S1、S2、S3，每个序列包含图像36张。图像经过预处理已将DICOM格式转换为易于计算机处理的PNG格式。

本实验的计算机配置为：cpu：Inter [email protected]，内存：4G，操作***Windows10，处理软件：Matlab2013a。

结果如图1至图5、表1、表2所示。图1为利用灰度分布信息对图像序列S1的均值图像进行ROI定位的效果图，添加了局部极小值搜索步骤，根据Matlab连通域提取函数获得的目标O′_k,m位置信息以及膝关节成像特点，可见此方法能准确的定位半月板并且无丢失的提取ROI。图2为利用灰度分布信息对图像序列的均值图像进行ROI定位的效果图，未添加局部极小值搜索步骤的结果对比图，图2a为未添加局部极小值搜索步骤时，利用灰度分布信息对图像序列S1的均值图像进行ROI定位的效果图，可见未添加局部极小值搜索步骤时ROI定位会由于成像因素的影响(膝关节在图像中的位置、膝关节弯曲程度等)，出现在图像两端，无法实现半月板准确定位，图2b、图2c分别为未添加局部极小值搜索步骤时，利用灰度分布信息对图像序列S2、S3的均值图像进行ROI定位的效果图，可见这种成像因素造成的影响并不是总会出现，但出现几率较高。图3a为序列S1中第5片层分割出的ROI图像序列示例图；图3b为序列S1中第7片层分割出的ROI图像序列示例图；图3c为序列S1中第9片层分割出的ROI图像序列示例图；图3d为序列S1中第11片层分割出的ROI图像序列示例图；图3e为序列S1中第13片层分割出的ROI图像序列示例图；可见ROI的选择可完整的包含半月板结构，并且使得面积尽可能减小便于运算、降低误提取概率。图4a为序列S1中第5片层半月板目标提取结果示例图；图4b为序列S1中第7片层半月板目标提取结果示例图；图4c为序列S1中第9片层半月板目标提取结果示例图；图4d为序列S1中第11片层半月板目标提取结果示例图；图4e为序列S1中第13片层半月板目标提取结果示例图；可见半月板结构被完整提取出来，并且已进行阈值化处理和空洞填充，此时半月板目标只有形状信息(失去位置、灰度信息)，从第11片层结果可见存在误提取现象(由于原目标较小不易展示，部分图像被放大)。图5a为序列S1中第5片层半月板信息提取结果示例图；图5b为序列S1中第7片层半月板信息提取结果示例图；图5c为序列S1中第9片层半月板信息提取结果示例图；图5d为序列S1中第11片层半月板信息提取结果示例图；图5e为序列S1中第13片层半月板信息提取结果示例图；可见半月板信息(包括形状、位置关系、灰度纹理信息)被完整地从原MRI中提取出来，图5a、图5b、图5c、图5d、图5e呈现的是整个ROI区域，从第11片层结果可见存在误提取现象(该部分为灰度相似，且满足形状约束的非半月板组织)。图6a为添加了基于医学先验知识结合位置关系和目标数量进行滤除的步骤后，序列S1中第11片层目标误提取现象消失对比图；图6b为添加了基于医学先验知识结合位置关系和目标数量进行滤除的步骤后，序列S1中第11片层目标误提取现象消失结果图。

表1半月板分割结果评估指标

表2不同形态半月板检测性能比较

ROI分割结果如图1所示，与未添加局部极小值搜索的方法结果对比如图2a-2c所示，可见实现了存在两端黑影干扰时的ROI准确定位。半月板目标提取结果如图3a-3e所示，与图4a-4e(未添加误提取目标滤除)对比，可见增加的步骤可降低误提取几率。图5a-5e显示了最终半月板分割的结果。利用专家对MRI序列中半月板目标的人工标注作为参考标准，计算出方法的性能指标：提取半月板目标的准确度、目标提取的召回率，如表1所示。由于方法中计算出的ROI面积以及阈值搜索范围小于已有算法的1/2，此部分运算效率为原有的二倍。实验中还分析了对于不同形状的半月板目标的提取性能，结果如表2所示。

附图说明

图1为添加了局部极小值搜索步骤时，利用灰度分布信息对图像序列S1的均值图像进行ROI定位的效果图；

图2a为未添加局部极小值搜索步骤时，利用灰度分布信息对图像序列S1的均值图像进行ROI定位的效果图；

图2b为未添加局部极小值搜索步骤时，利用灰度分布信息对图像序列S2的均值图像进行ROI定位的效果图；

图2c为未添加局部极小值搜索步骤时，利用灰度分布信息对图像序列S3的均值图像进行ROI定位的效果图；

图3a为序列S1中第5片层分割出的ROI图像序列示例图；

图3b为序列S1中第7片层分割出的ROI图像序列示例图；

图3c为序列S1中第9片层分割出的ROI图像序列示例图；

图3d为序列S1中第11片层分割出的ROI图像序列示例图；

图3e为序列S1中第13片层分割出的ROI图像序列示例图；

图4a为序列S1中第5片层半月板目标提取结果示例图；

图4b为序列S1中第7片层半月板目标提取结果示例图；

图4c为序列S1中第9片层半月板目标提取结果示例图；

图4d为序列S1中第11片层半月板目标提取结果示例图；

图4e为序列S1中第13片层半月板目标提取结果示例图；

图5a为序列S1中第5片层半月板信息提取结果示例图；

图5b为序列S1中第7片层半月板信息提取结果示例图；

图5c为序列S1中第9片层半月板信息提取结果示例图；

图5d为序列S1中第11片层半月板信息提取结果示例图；

图5e为序列S1中第13片层半月板信息提取结果示例图；

图6a为图4d经过滤除误提取目标处理后的对比图；

图6b为图5d经过滤除误提取目标处理后的对比图；

图7为子矩阵平均灰度曲线图，Mean intensity for submatrix为子矩阵的灰度平均值，Mean intensity为平均灰度，No.of submatrix为子矩阵的下标；

图8为归一化后的LSE的分布图，横坐标是n纵坐标是每个子矩阵计算的LSE；LSEbetween ROIs and mask为ROI和模板间的LSE，No.of ROIs为ROI的下标；

图9为本发明流程图；

图10为以均值化图像I_mean纵向中轴线为中心两侧各90像素宽度区域作为待处理区域，沿待处理区域纵向中轴线垂直方向分割为16行×180列像素大小的子矩阵示意图。

具体实施方式

具体实施方式一：结合图9说明本实施方式，本实施方式的一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，具体是按照以下步骤制备的：

可能包含半月板的片层可通过表3中信息初步选定，称为S_selected，所有片层称为S_all。

表3包含半月板的片层分布

半月板在膝关节MRI中出现的位置宽度，基本固定，通过经验可确定其宽度可选择在MRI中，中线左右宽度180像素的区域。

根据经验半月板高度一般在10-20像素之间。根据分析膝关节有可能出现的最大倾斜，通过计算可完全包含半月板的最小可接受ROI为_{60行×180列}(计算方法：假设以一种极端情况为例，两半月板上下错开50％厚度(10像素)，平均图像总高度达到30像素，再取极端情况，定位小片定位ROI在一个半月板远离另一个的一半处，则ROI一半宽度为30将完全包含该半月板5像素的厚度，另一半30像素，也足够包含另一边_15+15＝30)。

步骤一、ROI(感兴趣区域)定位，选定初始片层组S_selected；

所述，ROI为感兴趣区域；

步骤二、对选定的初始片层组S_selected进行优化，得到优化后片层集合S_optimized；

步骤三、对优化后片层集合S_optimized中的片层的ROI进行阈值自动选择与半月板目标提取；

步骤四、根据提取出的半月板目标在优化后片层集合S_optimized中进行半月板原始信息提取，得到包含分割出的半月板原始信息的ROI的片层T_k。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中ROI(感兴趣区域)定位，选定初始片层组S_selected；具体过程为：

(1)计算选定的初始片层的均值化图像，记为I_mean；

其中，所述选定为人为根据临床解剖学知识选定；选定的初始片层为512*512像素；

I_mean＝mean{I_k(i,j)|I_k∈S_selected} (1)

其中，I_k表示选定的初始片层组S_selected中的下标为k的片层；i为1-512像素横坐标，j为1-512像素纵坐标；

(2)以均值化图像I_mean纵向中轴线为中心两侧各90像素宽度区域作为待处理区域，沿待处理区域纵向中轴线垂直方向分割为16行×180列像素大小的子矩阵，如图10；

第n个子矩阵表示为P_n(i,j)，用于定位半月板垂直位置，子矩阵高度的选择接近根据经验的平均半月板厚度，同时与图像整体尺寸有关(512*512像素)；n取值为1-32；

(3)计算每个子矩阵的平均灰度值，如图7，并绘制出横坐标为n，纵坐标为灰度值的子矩阵平均灰度曲线图；

N_p(n)＝mean(P_n(i,j)),n∈(1,32) (2)

式中，N_p(n)为第n个子矩阵的平均灰度值；

(4)根据经验知识设置子矩阵的平均灰度值的局部极小值搜索范围为[30％，70％]，在范围为[30％，70％]内搜索平均灰度值最小值对应的子矩阵作为包含半月板小片；

$n^{*}=\arg \{\underset{n}{min}\[N_p\left(n\right)\]|n\∈(9,22)\}---\left(3\right)$

式中，n^*为平均灰度值最小的子矩阵的下标，n^*∈n；

将的纵向的高度/2+1行(的纵向的高度/2+1行为纵向的高度/2向下移动一行)作为ROI窗口的中心行，中心行向上选取30行，向下29行作为实际选定的ROI窗口，将实际选定的ROI窗口作用于全部片层获得ROI图像序列i取值为为1-180像素横坐标，j取值为为1-60像素纵坐标；

其中，为下标为k的片层的ROI，S_all为全部片层集合。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中对选定的初始片层组S_selected进行优化，得到优化后片层集合S_optimized；具体过程为：

(1)计算选定的初始片层组S_selected的ROI的均值作为模板φ(i,j)；

$\mathrm{\φ}(i,j)=\frac{1}{M}\underset{I_k\∈S_{selected}}{\Σ}{I}_k^{ROI}(i,j)---\left(4\right)$

其中，M为选定的初始片层组S_selected中的片层数；φ(i,j)为选定的初始片层组S_selected的ROI的均值；(2)利用模板φ(i,j)匹配所有选定片层I_k的ROI，计算最小均方误差LSE，并进行归一化处理，如图8；

$LSE\left(k\right)=\underset{i}{\Σ}\underset{j}{\Σ}{\left(I_k^{ROI}\right(i,j)-\mathrm{\φ}(i,j\left)\right)}^2,I_k\∈S_{selected}---\left(5\right)$

${\mathrm{LSE}}^{\′}\left(k\right)=\frac{LSE\left(k\right)}{\underset{k}{max}\left(LSE\right(k\left)\right)}---\left(6\right)$

将归一化处理后LSE′(k)>60％的片层视为不包含半月板片层，剔除这些片层，余下的片层作为优化后的片层集合S_optimized。

如果要选择出典型的包含类三角形半月板的MRI，则可将条件变严格，设置LSE′(k)>40％或者更小，则靠近两侧和中部的片层将被删除。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤三中对优化后片层集合S_optimized中的片层的ROI进行阈值自动选择与半月板目标提取；具体过程为：

(1)计算优化后片层集合S_optimized中的片层的ROI的均值图像的平均灰度，作为阈值选择的参考值；

$\mathrm{\η}=mean\left\{\underset{I_k\∈S_{optimized}}{\Σ}\underset{i}{\Σ}\underset{j}{\Σ}{I}_k^{ROI}(i,j)\right\}---\left(7\right)$

式中，η为优化后片层集合S_optimized中的片层的ROI的均值图像的平均灰度；

(2)选取阈值τ＝0.05η，对优化后片层集合S_optimized中的片层的ROI进行二值化处理；

$B_k(i,j)=\left\{\begin{array}{cc}255& I_k^{ROI}(i,j)\&GreaterEqual;\mathrm{\&tau;}\\ 0& I_k^{ROI}(i,j)<\mathrm{\&tau;}\end{array}\right.,I_k\&Element;S_{optimized}---\left(8\right)$

式中，B_k(i,j)为二值化处理后的图像；

(3)将二值化处理后的图像B_k(i,j)的ROI中所有满足八连通的独立目标利用Matlab连通域提取函数全部提取出来，进行滤波约束处理，记为目标O_k,m，k为片层编号，m为每片层中目标编号；此时目标包含位置坐标；

(4)首先进行过大或过小目标的滤除，将面积不满足表4中面积约束条件的目标滤除，初步减少噪声像素的干扰。

当目标O_k,m不满足面积为Area(O_k,m)>50像素，Area(O_k,m)<800像素时，滤除目标O_k,m；

表4目标形状约束条件

(5)针对滤除目标O_k,m后每一个剩余目标，计算剩余目标的轴比和边界盒比，轴比为主轴/副轴，边界盒比为边界长度/宽度，滤除不满足轴比和边界盒比的目标，得到满足半月板形状约束条件的目标O′_k,m，每个片层上ROI中目标数量记为下标t为阈值变化标记，1≤t≤19；形状约束条件见表4。

其中，轴比的目标为主轴/副轴>1，主轴/副轴<6，边界盒比目标为长度/宽度>1，长度/宽度<5；半月板形状约束条件为轴比目标为主轴/副轴>1，主轴/副轴<6，边界盒比目标为长度/宽度>1，长度/宽度<5；

$q_{k,t}^0=no.of\left(O_{k,m}^{\&prime;}\right)---\left(9\right)$

(6)根据Matlab连通域提取函数获得的目标O′_k,m位置信息以及膝关节成像特点，滤除误提取的目标，得到滤除误提取的目标后的提取结果O″_k，m(若出现剩余目标超过3个，则判定两侧目标为半月板，其余目标为关节软骨残余部分)；

$O_{k,m}^{\&prime;\&prime;}=\{\begin{array}{ccc}firstandlastof{O}_{k,m}^{\&prime;}& if& q_{k,t}^0>2\\ O_{k,m}^{\&prime;}& if& q_{k,t}^0\&le;2\end{array}---\left(10\right)$

式中，O″_k,m为滤除误提取的目标后的提取结果；

(7)记录每个片层上ROI中提取的半月板目标数量q_k,t；

q_k,t＝no.of(O″_k，m) (11)

以及阈值τ对应的ROI序列总共提取出的半月板目标数量；

$q_t=\underset{I_k\&Element;S_{optimized}}{\&Sigma;}{q}_{k,t}---\left(12\right)$

(8)改变阈值，重复步骤(3)—(7)，记录用每个阈值处理时提取出的半月板目标总数量q_t，阈值范围为τ＝0.05η～0.40η，每次增长0.05η，其中下标t＝1,2,......,8；

(9)选择使提取出的半月板目标总数量q_t最多的阈值作为自适应阈值选择的最优结果；

${\mathrm{\&tau;}}^{*}=\arg \left\{\underset{\mathrm{\&tau;}}{max}\left(q_t\right)\right\}---\left(13\right)$

以τ^*作为二值化处理的最优阈值，τ^*对应的提取结果O″_k,m即为最优半月板目标提取结果，此时提取半月板目标记为

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述τ^*取值为为0.10η～0.30η。通过实验我们发现对于我们的T1加权MRI数据，最优τ^*取值为0.10η～0.30η。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述步骤四中根据提取出的半月板目标在优化后片层集合S_optimized中进行半月板原始信息提取，得到包含分割出的半月板原始信息的ROI的片层T_k；具体过程为：

(1)根据半月板目标计算的凸包(填补阈值处理可能产生的目标内部空洞)作为模板。

(2)利用第二阶段步骤中Matlab连通域提取函数记录的每一个的位置坐标，获取半月板目标对应到原片层中的位置坐标。

(3)制作与ROI区域大小(60行*180列像素)相同的模板，其中对应于凸化处理后的半月板目标的区域为空，其他区域灰度设置为255(白色，便于展示与观察)；

T_k即为最终仅包含分割出的半月板原始信息的ROI的片层。

其它步骤及参数与具体实施方式一至五之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

本实施例一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法具体是按照以下步骤制备的：

实验中采用的三组实验数据是由哈尔滨医科大学附属医院提供，采集于三位年龄分布在30-50之间、于2014年在哈尔滨医科大学附属医院进行膝关节磁共振检查的患者，图像为T1加权(TR＝500ms，TE＝17ms)，片层间隔为4mm，像素尺寸为512*512。为尊重用户隐私，成像序列分别称为S1、S2、S3，每个序列包含图像36张。图像经过预处理已将DICOM格式转换为易于计算机处理的PNG格式。

本实验的计算机配置为：cpu：Inter [email protected]，内存：4G，操作***Windows10，处理软件：Matlab2013a。

结果如图1至图5、表1、表2所示。图1为利用灰度分布信息对图像序列S1的均值图像进行ROI定位的效果图，添加了局部极小值搜索步骤，根据Matlab连通域提取函数获得的目标O′_k,m位置信息以及膝关节成像特点，可见此方法能准确的定位半月板并且无丢失的提取ROI。图2为利用灰度分布信息对图像序列的均值图像进行ROI定位的效果图，未添加局部极小值搜索步骤的结果对比图，图2a为未添加局部极小值搜索步骤时，利用灰度分布信息对图像序列S1的均值图像进行ROI定位的效果图，可见未添加局部极小值搜索步骤时ROI定位会由于成像因素的影响(膝关节在图像中的位置、膝关节弯曲程度等)，出现在图像两端，无法实现半月板准确定位，图2b、图2c分别为未添加局部极小值搜索步骤时，利用灰度分布信息对图像序列S2、S3的均值图像进行ROI定位的效果图，可见这种成像因素造成的影响并不是总会出现，但出现几率较高。图3a为序列S1中第5片层分割出的ROI图像序列示例图；图3b为序列S1中第7片层分割出的ROI图像序列示例图；图3c为序列S1中第9片层分割出的ROI图像序列示例图；图3d为序列S1中第11片层分割出的ROI图像序列示例图；图3e为序列S1中第13片层分割出的ROI图像序列示例图；可见ROI的选择可完整的包含半月板结构，并且使得面积尽可能减小便于运算、降低误提取概率。图4a为序列S1中第5片层半月板目标提取结果示例图；图4b为序列S1中第7片层半月板目标提取结果示例图；图4c为序列S1中第9片层半月板目标提取结果示例图；图4d为序列S1中第11片层半月板目标提取结果示例图；图4e为序列S1中第13片层半月板目标提取结果示例图；可见半月板结构被完整提取出来，并且已进行阈值化处理和空洞填充，此时半月板目标只有形状信息(失去位置、灰度信息)，从第11片层结果可见存在误提取现象(由于原目标较小不易展示，部分图像被放大)。图5a为序列S1中第5片层半月板信息提取结果示例图；图5b为序列S1中第7片层半月板信息提取结果示例图；图5c为序列S1中第9片层半月板信息提取结果示例图；图5d为序列S1中第11片层半月板信息提取结果示例图；图5e为序列S1中第13片层半月板信息提取结果示例图；可见半月板信息(包括形状、位置关系、灰度纹理信息)被完整地从原MRI中提取出来，图5a、图5b、图5c、图5d、图5e呈现的是整个ROI区域，从第11片层结果可见存在误提取现象(该部分为灰度相似，且满足形状约束的非半月板组织)。图6a为添加了基于医学先验知识结合位置关系和目标数量进行滤除的步骤后，序列S1中第11片层目标误提取现象消失对比图；图6b为添加了基于医学先验知识结合位置关系和目标数量进行滤除的步骤后，序列S1中第11片层目标误提取现象消失结果图。

表1半月板分割结果评估指标

表2不同形态半月板检测性能比较

ROI分割结果如图1所示，与未添加局部极小值搜索的方法结果对比如图2所示，可见实现了存在两端黑影干扰时的ROI准确定位。半月板目标提取结果如图3a-3e所示，与图4a-4e(未添加误提取目标滤除)对比，可见增加的步骤可降低误提取几率。图5a-5e显示了最终半月板分割的结果。利用专家对MRI序列中半月板目标的人工标注作为参考标准，计算出方法的性能指标：提取半月板目标的准确度、目标提取的召回率，如表1所示。由于方法中计算出的ROI面积以及阈值搜索范围小于已有算法的1/2，此部分运算效率为原有的二倍。实验中还分析了对于不同形状的半月板目标的提取性能，结果如表2所示。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (6)

1.一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，其特征在于具体是按照以下步骤进行的：

步骤一、ROI定位，选定初始片层组S_selected；

所述，ROI为感兴趣区域；

步骤二、对选定的初始片层组S_selected进行优化，得到优化后片层集合S_optimized；

步骤三、对优化后片层集合S_optimized中的片层的ROI进行阈值自动选择与半月板目标提取；

步骤四、根据提取出的半月板目标在优化后片层集合S_optimized中进行半月板原始信息提取，得到包含分割出的半月板原始信息的ROI的片层T_k。

2.根据权利要求1所述一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，其特征在于：所述步骤一中ROI定位，选定初始片层组S_selected；具体过程为：

(1)计算选定的初始片层的均值化图像，记为I_mean；

选定的初始片层为512*512像素；

I_mean＝mean{I_k(i,j)|I_k∈S_selected} (1)

其中，I_k表示选定的初始片层组S_selected中的下标为k的片层；i为1-512像素横坐标，j为1-512像素纵坐标；

(2)以均值化图像I_mean纵向中轴线为中心两侧各90像素宽度区域作为待处理区域，沿待处理区域纵向中轴线垂直方向分割为16行×180列像素大小的子矩阵；

第n个子矩阵表示为P_n(i,j)；n取值为1-32；

(3)计算每个子矩阵的平均灰度值，

N_p(n)＝mean(P_n(i,j)),n∈(1,32) (2)

式中，N_p(n)为第n个子矩阵的平均灰度值；

(4)设置子矩阵的平均灰度值的局部极小值搜索范围为[30％，70％]，在范围为[30％，70％]内搜索平均灰度值最小值对应的子矩阵作为包含半月板小片；

$n^{*}=\arg \left\{\underset{n}{\min }\right[N_p\left(n\right)\left]\right|n\&Element;\left(\mathrm{9,22}\right)\}---\left(3\right)$

式中，n^*为平均灰度值最小的子矩阵的下标，n^*∈n；

将的纵向的高度/2+1行作为ROI窗口的中心行，中心行向上选取30行，向下29行作为实际选定的ROI窗口，将实际选定的ROI窗口作用于全部片层获得ROI图像序列I_k∈S_all，i取值为为1-180像素横坐标，j取值为为1-60像素纵坐标；

其中，为下标为k的片层的ROI，S_all为全部片层集合。

3.根据权利要求2所述一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，其特征在于：所述步骤二中对选定的初始片层组S_selected进行优化，得到优化后片层集合S_optimized；具体过程为：

(1)计算选定的初始片层组S_selected的ROI的均值作为模板φ(i,j)；

$\mathrm{\&phi;}(i,j)=\frac{1}{M}\underset{I_k\&Element;S_{selected}}{\&Sigma;}{I}_k^{ROI}(i,j)---\left(4\right)$

其中，

M为选定的初始片层组S_selected中的片层数；

φ(i,j)为选定的初始片层组S_selected的ROI的均值；

(2)利用模板φ(i,j)匹配所有选定片层I_k的ROI，计算最小均方误差LSE，并进行归一化处理；

$LSE\left(k\right)=\underset{i}{\&Sigma;}\underset{j}{\&Sigma;}{\left(I_k^{ROI}\right(i,j)-\mathrm{\&phi;}(i,j\left)\right)}^2,I_k\&Element;S_{selected}---\left(5\right)$

${\mathrm{LSE}}^{\&prime;}\left(k\right)=\frac{LSE\left(k\right)}{\underset{k}{max}\left(LSE\right(k\left)\right)}---\left(6\right)$

将归一化处理后LSE′(k)>60％的片层视为不包含半月板片层，剔除这些片层，余下的片层作为优化后的片层集合S_optimized。

4.根据权利要求3所述一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，其特征在于：步骤三中对优化后片层集合S_optimized中的片层的ROI进行阈值自动选择与半月板目标提取；具体过程为：

(1)计算优化后片层集合S_optimized中的片层的ROI的均值图像的平均灰度，作为阈值选择的参考值；

$\mathrm{\&eta;}=mean\left\{\underset{I_k\&Element;S_{optimized}}{\&Sigma;}\underset{i}{\&Sigma;}\underset{j}{\&Sigma;}{I}_k^{ROI}(i,j)\right\}---\left(7\right)$

式中，

η为优化后片层集合S_optimized中的片层的ROI的均值图像的平均灰度；

(2)选取阈值τ＝0.05η，对优化后片层集合S_optimized中的片层的ROI进行二值化处理；

$B_k(i,j)=\left\{\begin{array}{cc}255& I_k^{ROI}(i,j)\&GreaterEqual;\mathrm{\&tau;}\\ 0& I_k^{ROI}(i,j)<\mathrm{\&tau;}\end{array}\right.,I_k\&Element;S_{optimized}---\left(8\right)$

式中，

B_k(i,j)为二值化处理后的图像；

(3)将二值化处理后的图像B_k(i,j)的ROI中所有满足八连通的独立目标利用Matlab连通域提取函数全部提取出来，进行滤波约束处理，记为目标O_k,m，k为片层编号，m为每片层中目标编号；

(4)当目标O_k,m不满足面积为Area(O_k,m)>50像素，Area(O_k,m)<800像素时，滤除目标O_k,m；

(5)针对滤除目标O_k,m后每一个剩余目标，计算剩余目标的轴比和边界盒比，轴比为主轴/副轴，边界盒比为边界长度/宽度，滤除不满足轴比和边界盒比的目标，得到满足半月板形状约束条件的目标O′_k,m，每个片层上ROI中目标数量记为下标t为阈值变化标记，1≤t≤19；

其中，轴比的目标为主轴/副轴>1，主轴/副轴<6，边界盒比目标为长度/宽度>1，长度/宽度<5；半月板形状约束条件为轴比目标为主轴/副轴>1，主轴/副轴<6，边界盒比目标为长度/宽度>1，长度/宽度<5；

$q_{k,t}^0=no.of\left(O_{k,m}^{\&prime;}\right)---\left(9\right)$

(6)根据Matlab连通域提取函数获得的目标O′_k,m位置信息以及膝关节成像特点，滤除误提取的目标，得到滤除误提取的目标后的提取结果O″_k,m；

$O_{k,m}^{\&prime;\&prime;}=\left\{\begin{array}{ccc}firstandlastof{O}_{k,m}^{\&prime;}& if& q_{k,t}^0>2\\ O_{k,m}^{\&prime;}& if& q_{k,t}^0\&le;2\end{array}\right.---\left(10\right)$

式中，

O″_k,m为滤除误提取的目标后的提取结果；

(7)记录每个片层上ROI中提取的半月板目标数量_qk,t；

q_k,t＝no.of(O″_k,m) (11)

以及阈值τ对应的ROI序列总共提取出的半月板目标数量；

$q_t=\underset{I_k\&Element;S_{optimized}}{\&Sigma;}{q}_{k,t}---\left(12\right)$

(8)改变阈值，重复步骤(3)—(7)，记录用每个阈值处理时提取出的半月板目标总数量q_t，阈值范围为τ＝0.05η～0.40η，每次增长0.05η，其中下标t＝1,2,......,8；

(9)选择使提取出的半月板目标总数量q_t最多的阈值作为自适应阈值选择的最优结果；

${\mathrm{\&tau;}}^{*}=\arg \left\{\underset{\mathrm{\&tau;}}{max}\left(q_t\right)\right\}---\left(13\right)$

以τ^*作为二值化处理的最优阈值，τ^*对应的提取结果O″_k，m即为最优半月板目标提取结果，此时提取半月板目标记为

5.根据权利要求4所述一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，其特征在于：所述τ^*取值为为0.10η～0.30η。

6.根据权利要求5所述一种基于形状约束的膝关节磁共振图像序列半月板自动提取方法，其特征在于：所述步骤四中根据提取出的半月板目标在优化后片层集合S_optimized中进行半月板原始信息提取，得到包含分割出的半月板原始信息的ROI的片层T_k；具体过程为：

(1)根据半月板目标计算的凸包作为模板；

(2)利用Matlab连通域提取函数记录的每一个的位置坐标，获取半月板目标对应到原片层中的位置坐标；

(3)制作与ROI区域大小相同的模板，其中对应于凸化处理后的半月板目标的区域为空，其他区域灰度设置为255；

T_k即为最终仅包含分割出的半月板原始信息的ROI的片层。