一种医学图像中腹部器官分割方法

专利名称：一种医学图像中腹部器官分割方法
技术领域：
本发明属于医学图像处理及应用领域，涉及一种医学图像中腹部器官分割方法。
背景技术：
腹部器官的分割具有重要研究意义和临床价值。首先，从各种类型的医学图像中将器官分割出来是进行可视化的第一步。比可视化更重要的，确定感兴趣器官位置和区域在放疗手术计划中具有重要意义。只有掌握了器官的空间几何形状，才能制定出准确的放疗方案和放疗剂量。
Linda G. Shapiro等人1995年在《Pattern Recognition》杂志上提出利用器官形状的先验知识来动态调整灰度阈值进行分割，然后按照先易后难的顺序逐个识别感兴趣的器官。这种方法对每个器官采用的仍然是单一阈值分割，得到的区域往往是不连续的。John E.Koss等人1999年在《IEEE Transactions On Medical Imaging》杂志上提出选取图像的纹理特征，采用Hopfield网络进行聚类分割。这里分割得到的器官取决于Hopfield网络中指定的聚类数目。Chien-Cheng Lee等人2003年在《IEEE Transactions On Information Technology In Biomedicine》杂志提出采用上下文相关神经网络基础上的模糊识别规则对各区域进行识别和合并。这样一套模糊识别规则的建立，要求相当的经验和多次试验验证。Hyunjin Park等人2003年在《IEEETransactions On Medical Imaging》杂志提出则将关于腹部器官的解剖学形态图谱通过配准集成到贝叶斯分割框架中[5]。该方法需要手工设置多个控制点，影响了方法的自动化。

发明内容
本发明的目的是为临床应用提供一种在PC机上运行的腹部CT和核磁图像的器官分割算法，其特点是同时进行图像中腹部多器官自动化分割，实施简单，无需人工干预。
算法的主要步骤包括步骤(1)，在PC机中，采用基于归一化互信息的配准方法，将图谱向通过CT或核磁扫描得到的个体图像作整体配准，以消除图谱和目标图像间的整体差异；所述图谱是指对已有的医学图象经解剖学专家手工分割后对每个器官区域做好标记的图象，它可以为PC机提供人体解剖结构的参考信息；所述的归一化互信息nMI(S，T*A)用下式表示
nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>其中，S为目标图像，T*A代表对图谱A做空间变换T后的图像；H(S)代表目标图像的熵，H(T*A)代表变换后图谱的熵，H(S，T*A)代表目标图像和变换后图谱间的联合熵；nMI是英文的归一化互信息的缩写(normalized mutual information)，nMI(S，T*A)表示S与T*A之间的归一化互信息；空间变换T为包含9个参数的仿射变换T＝(p，q，r，u，v，w，φ，ω，θ)；其中，p，q，r分别为空间变换后x，y，z方向的位移；r，u，v分别为空间变换后x，y，z方向的比例缩放；φ，ω，θ分别为空间变换后绕x，y，z轴的转角；若变换前点坐标为X＝(x，y，z)，则变换后的坐标为X′＝X′，y′，z′)，x′y′z′1=100p010q001r0001×cosφsinφ00-sinφcosφ0000100001×00-sinω00100sinω0cosω00001]]>×10000cosθsinθ00-sinθcosθ00001×u0000v0000w00001×xyz1]]>至此，定义了目标图像S和经过空间变换T后的图谱A之间的归一化互信息；所述基于归一化互信息的配准方法依次按以下步骤实现步骤(1.1).随机确定初始的空间变换T0＝(T0，q0，r0，u0，v0，w0，φ0，ω0，θ0，)中九个参数值，并用由此得到的空间变换作用于图谱A，得到变换后的图谱T0*A；步骤(1.2).计算目标图像S和经过空间变换T后的图谱T*A间的联合归一化联合直方图h(l，k)h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}其中，S(x)＝l表示在目标图像中像素x对应的灰度值为l，T*A(x)＝k表示在经过空间变换的图谱T*A像素中x对应的灰度值为k；Num为像素x的数目；
步骤(1.3).计算目标图像S和经过空间变换T后的图谱T*A间的联合概率分布pS，T*A(l，K)pS,T*A(l,k)=h(l,k)Σl,kh(l,k).]]>步骤(1.4).分别计算图像S的边缘概率分布pS(l)和T*A的边缘概率分布pT*A(k)pS(l)=ΣkpS,T*A(l,k),]]>pT*A(k)=ΣlpS,T*A(l,k).]]>边缘概率分布是概率统计理论中的概念，对于两个随机变量x，y组成的随机向量(x，y)，其分量x的概率分布称为随机向量(x，y)的关于x的边缘分布；这里是把图像S和变换后的图谱T*A分别作为两个随机变量l和k，然后分别计算l和k的边缘概率分布；步骤(1.5).按以下公式分别计算目标图像的熵H(S)，变换后图谱的熵H(T*A)以及目标图像和变换后图谱间的联合熵H(S，T*A)；H(S)=-ΣlpS(l)logpS(l)]]>H(T*A)=-ΣkpT*A(k)logpT*A(k)]]>H(S,T*A)=-ΣkpS,T*A(l,k)logS,T*A(l,k)]]>步骤(1.6).按下式计算目标图像S和经过T变换后的图谱T*A间的归一化互信息nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>步骤(1.7).将归一化互信息nMI(S，T*A)带入Powell优化算法对空间变换T进行反复优化迭代，一直到收敛为止；Powell优化算法是一种传统的确定性优化方法，并且被认为是目前最有效的无需求导数的优方法；该算法的基本思想是，对于n维极值问题，首先沿着n个坐标方向求极小，经过n次之后得到n个共轭方向，然后沿n个共轭方向求极小，经过多次迭代后便可求得极小值；在这里，利用Powell优化算法对空间变换T的值进行优化，以求得到互信息nMI(S，T*A)的极小值；空间变换T含有九个分量，即(p，q，r，u，v，w，φ，ω，θ)，这九个分量分别对应Powell优化算法中的九个维度，即利用Powell优化算法求九维极值问题；当Powell优化算法搜索结束，即算法收敛后，得到的nMI(S，T*A)的极小值所对应的空间变换值即为最优空间变换T*，T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；则T**A为整体配准所要得到的配准图谱；步骤(2).依次按以下步骤对感兴趣的每一个器官进行单独配准，确定每个感兴趣的器官在目标图像S中的大致对应区域；具体分为以下几步步骤(2.1).按下式为每个感兴趣器官确定灰度范围区间RkRk＝[μk-λσk，μk+λσk]，其中，μk，σk为变换后图谱T*A的第k个感兴趣器官在目标图像S中所占据的区域内的灰度均值和标准方差；λ为一个常数，用来调整灰度均匀范围；试验发现，λ的取值范围为1.2～1.3；步骤(2.2).将步骤(1)中得到的最终的空间变换T*作对为每个器官单独进行配准的初始空间变换T′，即令T′＝T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；步骤(2.3).按下式计算变换后的器官区域的相似性测度Mk(T′)Mk(T′)=Nkin(T′)-Nkout(T′).]]>其中，下标k代表对第k个器官进行配准，Nkin(T′)代表采用空间变换T′时器官k在配准图谱中对应区域内包含的器官像素数目，Nkout(T′)代表该对应区域内包含的非器官像素数目；Nkin(T′)=N{x|x∈Ck(T′),S(x)∈Rk}Nkout(T′)=N{x|x∈Ck(T′),S(x)∉Rk}]]>其中，N{}表示计算集合中包含的像素数目，Ck(T′)表示器官k在空间变换T下对应的区域，S(x)表示像素x在CT或MRI图像中的灰度值；步骤(2.4).将变换后的器官区域的相似性测度Mk(T′)带入Powell优化算法对T′进行反复迭代，直到收敛为止；具体过程即利用Powell优化算法对空间变换T′的值进行优化，以求得到相似性测度Mk(T′)的极小值；空间变换T′含有九个分量，即(p′，q′，r′，u′，v′，w′，φ′，ω′，θ′)，这九个分量分别对应Powell优化算法中的九个维度，即利用Powell优化算法求九维极值问题；当Powell优化算法搜索结束，即算法收敛后，得到的Mk(T′)的极小值所对应的空间变换值即为最优空间变换T′*，T′*＝(p′*，q′*，r′*，u′*，v′*，w′*，ω′*，θ′*)；则Ck(T′*)为单独配准所要得到的器官区域；步骤(3).模糊连接分割；具体分为以下几步步骤(3.1).设Ck(T′*)为模糊连接对第k个器官进行分割的初始区域，按下式计算该区域中的灰度直方图Hk(l)Hk(l)＝N{x|x∈Ck，S(x)＝l}，l＝0，1……LHk(l)可理解为Ck(T′*)区域中灰度为l的像素x的个数；步骤(3.2).采用曲线拟合的方法对直方图进行平滑和去噪，得到拟合后的直方图 H^k(l)=Σj=1najexp[-(l-bjcj)2],]]>其中，n代表高斯分量的数目，即所要分割的器官的数目；例如我们对肝脏和脾脏两种器官进行分割，则n＝2；aj，bj，cj依次代表第j个高斯分布的峰值、中心和宽度；步骤(3.3).拟合后的直方图曲线可以看作器官在目标图像中的灰度概率密度函数；选择对应最大灰度概率的至少一个像素作为器官的初始种子点；步骤(3.4).考察每一个种子点相邻区域的灰度概率之和；相邻区域定义为与作为中心像素的种子点距离小于一个像素单位的距离的像素集合；所述灰度概率为所研究器官区域内具有该灰度的像素数目与区域内像素总数的比值；具有最大区域灰度概率总和的点将被最终确定为种子点；步骤(3.5).以由步骤(3.2)得到的直方图 所对应的器官区域灰度分布的概率密度函数和最终种子点的位置为初始参数，用下式所述的用模糊连接算法计算整幅图像中每一像素点与器官最终种子点的模糊连接强度cμ(p,q)=maxρ(p,q)∈P(p,q)[minz∈ρ(p,q)Hk^(l(z))]]]>其中，p表示最终的种子点，q表示整幅图像中任意一点，cμ(p，q)表示p与q之间的模糊连接强度；ρ(p，q)表示p与q之间的任意一条连通路径，P(p，q)表示p与q之间所有连通路径的集合，z∈ρ(p，q)表示Z是路径ρ(p，q)上的任意一像素点，l(z)为像素Z在目标图像中的灰度值， 表示灰度值l(z)在直方图 中所对应的概率密度值；模糊连接强度表征的该像素属于器官区域的隶属度，将每个像素点的隶属度作为灰度所得到的图像称为模糊连接图像；步骤(3.6).为每个器官的模糊连接图像指定一个最佳阈值；所谓阈值，就是用来将属于目标器官的像素和背景像素区分开的一个数值，规定灰度大于或等于此阈值的像素为属于目标器官的像素，灰度小于此阈值的像素为背景像素；所谓最佳阈值，就是由此阈值得到的器官像素所集合成的区域最符合真实的器官区域；最佳阈值是将模糊连接图像转化为最终的分割结果的重要参量；确定最佳阈值具体步骤如下首先将阈值设定为1，然后逐步降低阈值；用于每一次降低得到的阈值来分割模糊连接图像，得到器官的分割区域；观察分割区域的形状随阈值的改变，如果在某一次阈值降低时，形状剧烈改变，则停止降低阈值，将上一阈值定为最佳阈值；关于分割区域的形状改变的衡量办法和最佳阈值的确定方法，进一步详细解释如下a).分割区域的形状改变的衡量的办法采用面积相对变化值和紧密度相对变化值对来衡量分割区域的形状改变；紧密度(Compactness)的定义为 面积相对变化值和紧密度相对变化值定义如下对于每次阈值降低，分别计算降低前后的器官面积和紧密度；令降低前的器官面积和紧密度为原面积和原紧密度，降低后的器官面积和紧密度为新面积和新紧密度，则面积和紧密度的相对变化定义为 b).用面积相对变化值和紧密度相对变化值确定最优阈值的方法在阈值降低过程中，如果在某一次降低时，发现面积和紧密度的相对变化值大于预先设定的参考值ρ，则认为形状发生了剧烈改变，从而停止阈值搜索，并将降低前的阈值定为最佳阈值；试验表明，参考值ρ取值0.15～0.2能够取得较准确的最佳阈值；步骤(3.7)、按照下式用步骤(3.6)得到的最佳阈值对模糊连接图像进行分割，得到器官的二值分割图像设第k个器官的最佳阈值为Tk，则该器官的二值分割结果为 其中，x代表任一像素，Fk(x)为第k个器官的模糊连接图像在像素x处的灰度值，Bk(x)为分割后的二值图像在像素x处的灰度值；步骤(4).器官形状修正，以解决因为人体腹腔内相互紧贴的两个器官因为灰度相近而造成模糊连接分割可能产生的错误分割结果；医学图像中，经常出现两个在人体腹腔中相互紧贴的器官在图像中的灰度也相似的现象；对于这样的两个器官，医学图像中在其表面相贴处可能没有明显的边界，这会使模糊连接算法认为这两个器官是长在一起的，从而错误地分割成同一个器官；本研究提出以器官形状稳定为依据来检测错误分割，并采用距离变换结合分水岭分割算法的对错误分割结果进行形状修正；具体流程如下步骤(4.1).对于多层医学图像，按以下步骤从上到下，采用逐层分析的方法检测存在错误分割的层面；步骤(4.1.1).
设定当前层为算法正在分析的层面，当前层的序号为i下方相邻层为与当前层相邻的下方层面，下方相邻层的序号为i-1；Oki为当前层中的器官分割区域；
Oki-1为下方相邻层中器官分割区域；在当前层中，令器官分割区域内的各点到区域边界的距离为负值，对于整幅图像中所有不属于该区域的点，即所有该区域以外的点，令其到区域边界的距离为正值；步骤(4.1.2).对当前层i中的器官分割区域Oki进行距离变换，即对整幅图像中每点按下式计算它们到区域边界的最近距离Dki 其中，y代表了区域边界的像素点中与x距离最近一点；Bki(x)为步骤(3.7)中得到的最终分割结果Bk*(x)在当前层i中的像素值；将距离值dki(x)作为灰度值赋给每一像素点，就得到了距离图像Dki；步骤(4.1.3).提取下方相邻层中区域Oki-1的边界上各点像素的位置，则在Dki中的相同位置处的各点的Dki(x)值可认为是Oki-1边界上的像素的距离值；得到所有Oki-1的边界像素的距离值，计算这些距离值的标准方差，以确定器官在第i层和第i-1层分割结果间的相似程度；如果标准方差大于预先设定的参考值，则认为检测到了分割错误，下一层即为存在错误分割的层面；根据经验，这个参考值一般设定在5个像素距离；按照上述方法，从上到下检测所有层面，得到所有存在错误分割的层面；步骤(4.2).按下式对存在错误分割的各层面的二值图像进行距离变换；设定对于器官区域内所有像素点，令其距离值为负值；对于整幅图像中区域外所有不属于该区域的区域以外的点，令其距离值为负的极小值，即负无穷-∞；则距离图Dki′(x) 其中y代表了与x距离最近的边界点；Bki即存在错误分割的层面的二值图像；步骤(4.3).按以下方法对步骤(4.2)得到的距离图像Dki′应用分水岭算法；步骤(4.3.1).对整幅图像设定一个灰度阈值T；令T的初值为图像Dki′中除了负无穷以外的最小值，即
T0=min({Dki(x)|Dki(x)≠-∞}),]]>其中，x代表图像中任一个像素点的坐标；步骤(4.3.2).以d为步长逐渐增大阈值T，即Tn＝T0+n·d，其中，Tn表示第n次增加步长时阈值T的取值；n取值范围1≤n≤K，其中K=(max(Dki(x))-T0)/d,]]>max(Dki(x))表示图像Dki′的最大值；d取值固定，根据经验取d＝3；步骤(4.3.3).利用第n次增加步长得到的阈值Tn按下式将图像Dki′分割成二值图像Wn 其中，x代表图像中任一个像素点的坐标；Wn(x)＝1表示淹没在水面以下的部分，Wn(x)＝0表示露在水面以上的陆地部分；对于前一次(即第n-1次)增加步长得到阈值Tn-1，按下式将图像Dki′分割成二值图像Wn-1 其中，x代表图像中任一个像素点的坐标；步骤(4.3.4).在Tn的增长过程中， 每当Tn大于一个局部极大值，就说明图像Wn-1内有两个连通区域在Wn中合并成了一个连通域，则将Wn-1和Wn两幅图像进行异或操作，得到分水岭图像Sn 其中，x代表图像中任一个像素点的坐标，表示异或操作，Sn(x)＝1表示像素x为分水岭所在位置；步骤(4.3.5).对n取值范围1≤n≤K中的每一个n值，均会得到一个相应的分水岭图像Sn；设所有的n值所对应的分水岭图像为S1，S2，…， 则分水岭变换的结果图像S为所有分水岭图像的集合
S=S1·S2·...·Smax(Dki(x));]]>图像S中所有像素值为1的像素就表示图像Dki′中的分水岭所在的位置；由这些分水岭即把图像Dki′分割成若干个子区域；步骤(4.4).对分水岭算法得到的图像Dki′中的每个子区域，分别计算其与下方相邻层面器官分割区域Oki-1的面积重合比率 把重合比率低于50％的区域视为错误分割区域，重新标记为背景区域，从而将其在距离图像Dki′中除去。
实验效果分析为了验证本分割算法的可靠性，实验选择了5套来自于国内医院的临床CT腹部扫描图像。这5套数据包括男性和女性，以及不同年龄的对象。从分割结果看出，试验数据中的感兴趣器官均得到了有效分割。配准后图谱给出了大致分布区域，经过进一步的分割和修正，最终分割结果和手工分割结果相比，二者基本一致。


图1为发明采用的VIP-Man图谱中的一个层面。
图2为CT图像的一个层面。
图3为计算机程序的主体流程图。
图4为将图谱向CT图像整体配准后，提取图谱中感兴趣器官的轮廓叠加显示在CT图像上。感兴趣器官包括肝脏、肾脏和脾脏。
图5为器官配准后，图谱中感兴趣器官的轮廓和CT图像中器官位置对比。
图6为各个器官的灰度直方图和曲线拟合示意图。
图7为肝脏的模糊连接分割过程。从左至右依次为，CT图像中肝脏局部，肝脏的模糊连接图像，对模糊连接图像进行阈值分割后的肝脏分割结果。
图8为分水岭算法示意图。
图9为对肾脏进行形状修正过程。从左至右依次为，CT图像中肾脏局部，肾脏的模糊连接分割结果，肾脏经过分水岭算法被分割为多个区域和白色轮廓表示的相邻层面肾脏区域对比，消除重合比率过低的区域后得到的肾脏修正后分割结果。
图10对CT图像的整体分割结果。
具体实施例方式
本发明提出了一种基于图谱匹配的腹部器官分割方法，实现医学图像中腹部多器官的自动化分割。发明采用由美国虚拟人CT数据经过手工分割建立的图谱，对腹部器官进行分割。关于该图谱建立方法的详细说明可以参见Xu等人2000年在《Health Physics》杂志发表的论文《VIP-manAn image-based whole-body adult male model constructed from color photographs ofthe visible human project for multi-particle Monte Carlo calculations》。该图谱的示例可见图1。通过配准工作，得到的是一个和目标图像大致对齐的图谱，目标图像即需要做器官分割的临床图像，其示例见图2。配准的目的就是将图谱通过形变和目标图像大致对齐，这种对齐集中于感兴趣器官的对齐，而不是图像所有区域的对齐，这样就大大减轻了配准工作的压力和困难程度。从而允许在配准中选择较简单的空间变换，只包含较少的自由度。优化过程也相应的提高了速度，增强了稳定性。在分割算法上，选择模糊连接分割算法，它具有较强的稳定性，能够适应复杂图像的分割。模糊连接分割面临的一个主要问题是各项必要参数的指定，包括初始的种子点，感兴趣器官的灰度分布特征等。利用配准后图谱，本发明提出了这些参数的自动化指定方法，避免了常用的人工或经验指定。已发现，由于相邻器官间的灰度重叠或重合现象，模糊连接分割的结果可能包含着错误的区域。为了消除这些错误区域，提高分割结果的准确率，本研究基于器官在上下两张相邻层面的CT图像中具有形状一致性的原理，提出了器官形状修正方法，对分割结果做进一步的完善。
本方法的计算机程序的主体流程图可参加图3。
本发明方法通过下述步骤实现1.整体配准整体配准的目标是消除图谱和目标图像间的整体差异，包括成像位置以及身体形态(高矮胖瘦等)的差异。整体配准，寻找的是一种空间变换，使图谱经过这一空间变换后，能够在整体上与目标图像大致对齐。整体配准中选择相似变换(包含相对于x、y、z轴3个坐标轴的平移、旋转和比例变换，共9个参数)，以归一化的互信息作为目标函数，采用Powell算法进行优化。
图谱和目标图像间归一化互信息(Normalized Mutual Information，以下缩写为nMI)定义为nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>其中A代表图谱，T*A代表对图谱做空间变换T后的图像，H(S)代表目标图像的熵，H(T*A)代表变换后图谱的熵，H(S，T*A)代表目标图像和变换后图谱间的联合熵。
在空间变换上，选择包含9个参数的仿射变换T＝(p，q，r，u，v，w，φ，ω，θ)。其中，p，q，r分别为空间变换后x，y，z方向的位移；u，v，w分别为空间变换后x，y，z方向的比例缩放；φ，ω，θ分别为空间变换后饶x，y，z轴的转角；令变换前点坐标为X＝(x，y，z)，变换后的坐标为X′＝(x′，y′，z′)，则x′y′z′1=100p010q001r0001×cosφsinφ00-sinφcosφ0000100001×00-sinω00100sinω0cosω00001]]>×10000cosθsinθ00-sinθcosθ00001×u0000v0000w00001×xyz1]]>H(S)=-ΣlpS(l)logpS(l)]]>H(T*A)=-ΣkpT*A(k)logpT*A(k)]]>H(S,T*A)=-ΣkpS,T*A(l,k)logS,T*A(l,k)]]>这里，pS，T*A(l，k)是图像S和T*A间的联合概率分布，可以用归一化的联合直方图表示pS,T*A(l,k)=h(l,k)Σl,kh(l,k).]]>ps(l)，pT*A(k)是图像S和T*A间的边缘概率分布，由以下公式计算pS(l)=ΣkpS,T*A(l,k),]]>pT*A(k)=ΣlpS,T*A(l,k).]]>边缘概率分布是概率统计理论中的概念，对于两个随机变量x，y组成的随机向量(x，y)，其分量x的概率分布称为随机向量(x，y)的关于x的边缘分布；这里是把图像S和变换后的图谱T*A分别作为两个随机变量l和k，然后分别计算l和k的边缘概率分布；归一化联合直方图h(l，k)定义为
h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}，即，像素在图像S中灰度为l，并且对应像素在图像T*A中为k的像素数目统计。
至此，定义了目标图像S和经过空间变换T后的图谱A之间的互信息。
Powell优化算法是一种传统的确定性优化方法，并且被认为是目前最有效的无需求导数的优方法。该算法的基本思想是，对于n维极值问题，首先沿着n个坐标方向求极小，经过n次之后得到n个共轭方向，然后沿n个共轭方向求极小，经过多次迭代后便可求得极小值。
在这里，利用Powell优化算法对空间变换T的值进行优化，以求得到互信息nMI(S，T*A)的极小值。空间变换T含有九个分量，即(p，q，r，u，v，w，φ，ω，θ)，这九个分量分别对应Powell优化算法中的九个维度，即利用Powell优化算法求九维极值问题。当Powell优化算法搜索结束，即算法收敛后，得到的nMI(S，T*A)的极小值所对应的空间变换值即为最优空间变换T*，T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；则T**A为整体配准所要得到的配准图谱；图4展示了将图谱向CT图像整体配准后的效果。我们将经过最优空间T*变换后的图谱中的感兴趣器官的轮廓提取出来并叠加显示在CT图像上，以示效果。
2.器官配准器官配准的基本思想是进一步消除整体配准中未能解决的对应器官间的位置差异。基本方法是对图谱中每个感兴趣器官分别地向个体图像上配准，同时还要保证各个器官间不发生位置冲突。
对各个器官的单独配准仍然选择相似变换，采用Powell优化算法。器官配准定义一个全新的目标函数，图谱中的每一个感兴趣器官在目标图像中都存在一个对应区域，器官配准将努力使得该对应区域内包含尽可能多的器官像素，和尽可能少的非器官像素。由于这里还无法确定每个像素是否属于某个器官，因此，采用一种推测标准。设共有n个感兴趣器官，则对其中第k个感兴趣器官，其灰度区间推测为，Rk＝[μk-λσk，μk+λσk]，其中，μk，σk为图谱的第k个感兴趣器官经配准后所占据的区域内的灰度均值和标准方差。
其中λ为一个常数，用来调整灰度均匀范围。试验发现，λ的取值范围为1.2～1.3。如像素灰度在区间Rk内，则认为其属于器官k像素；否则，认为其属于非器官像素。配准的目标是通过最优化算法，对每一个器官分别确定一个新的相似变换；这里采用的最优化算法为Powell算法。新的变换使得器官对应区域内器官像素尽可能多，非器官像素尽可能少。定义器官配准的相似性测度Mk(T)如下Mk(T)=Nkin(T)-Nkout(T).]]>其中，下标k代表对第k个器官进行配准，Nkin(T)代表采用空间变换T时器官k的对应区域内包含的器官像素数目，Nkout(T)代表对应区域内包含的非器官像素数目。
Nkin(T)=N{x|x∈Ck(T),S(x)∈Rk}Nkout(T)=N{x|x∈Ck(T),S(x)∉Rk}]]>其中，N{}表示计算集合中包含的像素数目，Ck(T)表示器官k在空间变换T下对应的初始区域，S(x)表示像素在CT或MRI图像中的灰度值。
Powell优化算法的具体优化过程为将变换后的器官区域的相似性测度Mk(T′)带入Powell优化算法，利用Powell优化算法对空间变换T′的值进行优化，以求得到相似性测度Mk(T′)的极小值；空间变换T′含有九个分量，即(p′，q′，r′，u′，v′，w′，φ′，ω′，θ′)，这九个分量分别对应Powell优化算法中的九个维度，即利用Powell优化算法求九维极值问题；当Powell优化算法搜索结束，即算法收敛后，得到的Mk(T′)的极小值所对应的空间变换值即为最优空间变换T′*，T′*＝(p′*，q′*，r′*，u′*，v′*，w′*，φ′*，ω′*，θ′*)；则Ck(T′*)为单独配准所要得到的器官区域；图5展示了经过器官配准后，图谱中的感兴趣器官的轮廓和CT图像中器官位置对比。
3.模糊连接分割模糊连接算法由J.Udupa提出，它通过定义图像中每一像素点与器官种子点的模糊连接强度来得到每个像素属于该器官的隶属度。将每个像素点的隶属度作为灰度的图像即模糊连接图像。为了将模糊连接方法应用于器官分割，首先需要指定每个感兴趣器官的灰度分布特征和种子点位置(器官的种子点是指模糊连接算法分割开始分割时的初始像素点)。这些参数确定之后，模糊连接算法会分别计算每个感兴趣器官的模糊连接图像；模糊连接图像是灰度图像，其中每个像素点的灰度值代表了该点与种子点的模糊连接程度。模糊连接图像生成后，还需要进一步指定一个阈值，将其转化为对器官的二值分割图像，作为最终的分割结果。
配准后的图谱给出了器官在目标图像中的大致位置和范围，称其为器官的初始区域。统计初始区域中像素的灰度直方图。以第k个器官为例，灰度直方图的计算公式为，Hk(l)＝N{x | x∈Ck，S(x)＝l}， l＝0，1……L，即在器官k的初始区域中灰度为l的像素的数目。
器官灰度直方图反映了像素灰度值与像素属于器官的可能性大小之间的关系。由于初始区域中包含的不全部是器官像素，还有非器官像素。因此采用曲线拟合的方法对直方图进行平滑和去噪。曲线拟合将直方图视为若干个高斯函数的叠加，采用最小方差进行拟合。
器官拟合公式为H^k(l)=Σj=1najexp[-(l-bjcj)2],]]>其中，n代表高斯分量的数目，即所要分割的器官的数目；例如我们对肝脏和脾脏两种器官进行分割，则n＝2；aj，bj，cj依次代表第j个高斯分布的峰值、中心和宽度。在高斯分量数目和拟合精度之间，需要进行平衡，用尽可能少的高斯分量取得尽可能高的拟合精度。这里，采用调整的残值平方(Adjusted R-square，简称AR)来寻找最优的高斯分量数目。
AR=1-SSE*LSST*(L-n*3),]]>其中SSE=Σl=0L(Hk(l)-H^k(l))2]]>SST=Σl=0L(Hk(l)-H‾k)2,]]>H‾ki=Σl=0LHk(l)/(L+1).]]>具体的高斯拟合过程为(1)令n＝0；(2)令n＝n+1；(3)采用最小SSE方法，拟合式 (4)计算AR；(5)当AR＜AR并且n＜5返回执行步骤(2)；
在上述步骤中，TAR是依据经验设定的阈值，用于控制平衡拟合精度和高斯分量数目。试验表明，TAR＝0.95～0.97是一个比较合理的阈值设定。
图6展示了各个器官的灰度直方图和曲线拟合示意图。拟合后的直方图曲线可以看作器官在目标图像中的灰度概率密度函数。选择对应最大灰度概率的至少一个像素作为器官的初始种子点。对应最大灰度概率的像素数目可能大于1，再进一步考察像素相邻区域的灰度概率之和。相邻区域一般定义为与中心像素距离小于一个像素单位的像素集合。所述区域灰度概率为所研究器官区域内具有该灰度的像素数目与区域内像素总数的比值；具有最大区域灰度概率总和的点将被最终确定为种子点。同时规定，种子点必须位于器官的初始区域内部。
确定了器官的灰度分布特征和初始种子点，可以计算出每一个器官各自的模糊连接图像。以上面得到的直方图 所对应的器官区域灰度分布的概率密度函数和最终种子点的位置为初始参数，用下式所述的用模糊连接算法计算整幅图像中每一像素点与器官最终种子点的模糊连接强度cμ(p,q)=maxρ(p,q)∈P(p,q)[minz∈ρ(p,q)H^k(l(z))]]]>其中p表示最终的种子点，q表示整幅图像中任意一点，cμ(p，q)表示p与q之间的模糊连接强度；ρ(p，q)表示p与q之间的任意一条连通路径，P(p，q)表示p与q之间所有连通路径的集合，z∈ρ(p，q)表示Z是路径ρ(p，q)上的任意一像素点，l(z)为像素Z在目标图像中的灰度值， 表示灰度值l(z)在直方图 中所对应的概率密度值；模糊连接强度表征的该像素属于器官区域的隶属度，将每个像素点的隶属度作为灰度所得到的图像称为模糊连接图像。而要将器官真正分割出来，还需要指定一个模糊连接强度阈值，对模糊连接图像进行二值分割。所谓阈值，就是用来将目标像素和背景像素区分开的一个数值；在模糊连接图像中，我们认为灰度值大于或等于阈值的像素为器官像素，灰度值小于阈值像素为背景像素。所谓最佳阈值，就是由此阈值得到的器官像素所集合成的区域最符合真实的器官区域；最佳阈值是将模糊连接图像转化为最终的分割结果的重要参量。
为了找到最佳的模糊连接分割阈值，提出如下的搜索方法。首先将阈值设定为1，然后逐步降低阈值；用于每一次降低得到的阈值来分割模糊连接图像，得到器官的分割区域；分割方法如下设第k个器官的阈值为Tk，则该器官的阈值分割结果为
其中x代表任一像素，Fk(x)为第k个器官的模糊连接图像在像素x处的灰度值。
对于每一次阈值降低，观察分割结果随阈值的形状改变。如果在某一次阈值降低时，形状剧烈改变，则停止阈值搜索，将上一阈值定为最佳阈值。这里对形状改变的衡量使用了两个参量，即面积相对变化值和紧密度相对变化值。紧密度(Compactness)的定义为 对于每次阈值降低，分别计算降低前后的器官面积和紧密度；令降低前的器官面积和紧密度称为原面积和原紧密度，降低后的器官面积和紧密度称为新面积和新紧密度，则面积和紧密度的相对变化定义为 如果在某一次阈值降低时，发现面积和紧密度的相对变化值大于预先设定的参考值ρ，则认为形状发生了剧烈改变，从而停止阈值搜索，并将降低前的阈值定为最佳阈值。试验表明，参考值ρ取值0.15～0.2能够取得较准确的最佳阈值。
接下来，用上面确定的最佳阈值将模糊连接图像分割为二值图像，作为最终的输出。具体方法如下设第k个器官的最佳阈值为Tk*，则该器官的二值分割结果为 其中x代表任一像素，Fk(x)为第k个器官的模糊连接图像在像素x处的灰度值，Bk*(x)为分割后的二值图像在像素x处的灰度值。
图7以肝脏为例展示了模糊连接分割过程。从左至右依次为，CT图像中肝脏局部，肝脏的模糊连接图像，对模糊连接图像进行阈值分割后的肝脏分割结果。
4.器官形状修正医学图像中，经常出现两个在人体腹腔中相互紧贴的器官在图像中的灰度也相似的现象；对于这样的两个器官，医学图像中在其表面相贴处可能没有明显的边界，这会使模糊连接算法认为这两个器官是长在一起的，从而错误地分割成同一个器官；本研究提出以器官形状稳定为依据来检测错误分割，并采用距离变换结合分水岭分割算法的对错误分割结果进行形状修正；具体流程如下
1).对于多层医学图像，采用逐层分析的方法检测存在错误分割的层面；检测顺序从最上方层面到最下方层面；对于每次层面移动，称算法正在分析的层面为当前层，与当前层相邻的下方层面称为下方相邻层；设当前层的序号为i，则下方相邻层的序号为i-1；在当前层中，设Oki代表器官分割区域；对Oki进行距离变换，即计算区域内每个离散点到区域边界的距离；对于区域内所有点，令其距离值为负值；对于区域外所有点(范围包括整幅图像中所有不属于该区域的点)，令其距离值为正值，这样得到距离图像为Dki， 其中x代表图像中任意一处的离散像素点，y代表了区域边界上与x距离最近的点；Bki(x)为模糊连接最终分割结果Bk*(x)在第当前层的像素值。
在下方相邻层中，假定Oki-1代表器官分割区域；提取Oki-1的边界像素的位置，则在Dki中的相同位置处的Dki(x)值可认为是Oki-1边界上的像素与Oki边界上的对应像素的距离值；当得到所有Oki-1的边界像素的这种距离值后，计算这些距离值的标准方差，该标准方差反映了器官在第i层和第i-1层分割结果间的相似程度。如果器官从一层到另一层只是均匀的膨胀或收缩，则标准方差取值应小于预定义参考值。因此，可以采用器官在相邻层面间边界距离标准方差的大小来衡量是否出现了错误的异常突出区域。如果标准方差大于参考值，则认为检测到了分割错误。根据试验经验，这个参考值一般设定在5个像素距离。
按照上述方法，从上到下检测所有层面，得到所有存在错误分割的层面。
2).对存在错误分割的层面的二值图像进行距离变换；这里的距离变换与上面有所不同；对于区域内所有点，令其距离值为负值；对于区域外所有点(范围包括整幅图像中所有不属于该区域的点)，令其距离值负的极小值；这样得到距离图像Dki′为 其中x代表图像中任意一处的离散像素点，y代表了区域边界上与x距离最近的点；Bki即存在错误分割的层面的二值图像；3).对图像Dki′应用分水岭算法。
借助涨水法来介绍分水岭技术。分水岭(Watershed)算法把图像中像素的灰度值看作地形图上的海拔高度，则图像中灰度高的部分对应山峰，灰度低的部分对应山谷。利用图8(a)来进行讨论，为了简便，仅画出图像的一维剖面图。横坐标代表图像中的像素点，纵坐标代表像素的灰度值。假设有水从谷底涌出并水位逐渐升高。如果从两个谷底涌出的水水位高过了其间的山峰，这些水就会汇合。如果要阻止水的汇合，则需要在山峰上修建水坝，且水坝的高度随水位的上升而增高。这个过程随着全部山峰都被水淹没而结束。修建的所有水坝把整幅地形图分割成很多区域，这些水坝就构成了这片土地的分水岭。被分水岭分割出来的每一个区域对应于图像中的一个子区域，因此分水岭算法将距离图像Dki′分成多个子区域。用分水岭算法将距离图像Dki′分成多个子区域的具体过程如下仍以图8(a)而言，图中出现的三个山峰表示从图像Dki′中任意取出的一部分，山峰表示灰度高的像素值，山谷表示灰度低的像素值。设灰度阈值T代表水面的高度；令T的初值为图像Dki′中除负无穷以外的最小值，即T0=min({Dki(x)|Dki(x)≠-∞}),]]>其中，x代表图像中任一个像素点的坐标。
以d为步长逐渐增大阈值T，即Tn＝T0+n·d，其中，Tn表示第n次增加步长时阈值T的取值；n取值范围1≤n≤K；K=(max(Dki(x))-T0)/d,]]>max(Dki(x))表示图像Dki′的最大值；d取值固定，根据经验取d＝3。
利用第n次增加步长得到的阈值Tn(如图8(c)所示)，可以按下式将图像Dki′分割成二值图像Wn 其中，x代表图像中任一个像素点的坐标。Wn(x)＝1表示淹没在水面以下的部分，Wn(x)＝0表示露在水面以上的陆地部分。
对于前一次(即第n-1次)增加步长得到阈值Tn-1(如图8(b)所示)，可以按下式将图像Dki′分割成二值图像Wn-1
其中，x代表图像中任一个像素点的坐标。
在Tn的增长过程中，每当Tn大于一个局部极大值，就说明图像Wn-1内的两个连通区域(如图8(b)中的连通域Cj和Cj+1)在Wn中合并成了一个连通域(如图8(c)中的Cj，j+1)，则将Wn-1和Wn两幅图像进行异或操作，得到分水岭图像Sn 其中，x代表图像中任一个像素点的坐标，表示异或操作，Sn(x)＝1表示像素x为分水岭所在位置，如图8(d)所示。
因为n取值范围1≤n≤ K，对于这个范围的每一个n值，均会得到一个相应的分水岭图像Sn；则所有的n值所对应的分水岭图像为S1，S2，…， 分水岭变换的结果图像S为所有分水岭图像的集合S=S1·S2·...·Smax(Dki(x)).]]>则图像S中所有像素值为1的像素就表示图像Dki′中的分水岭所在的位置；由这些分水岭即把图像Dki′分割成若干个子区域。
4).对分水岭算法得到的每个子区域，分别计算其与下方相邻层面器官分割区域Oki-1的面积重合比率；面积重合比率定义为 基于器官在相邻层面中形状一致的原则，认为重合比率低于50％的区域为应为错误分割区域，应重新标记为背景区域，，从而将其在距离图像Dki′中除去。
图9以肾脏为例说明了形状修正过程。从左至右依次为，CT图像中肾脏局部，肾脏的模糊连接分割结果，肾脏经过分水岭算法被分割为多个区域和白色轮廓表示的下方相邻层面肾脏区域对比，消除重合比率过低的区域后得到的肾脏修正后分割结果。
至此，就全部完成了对腹部医学图像的分割，最终的分割结果可见图10。
权利要求
1.一种医学图像中腹部器官分割的方法，具体特征在于，该方法是在PC机中依次按以下步骤实现的算法的主要步骤包括步骤(1)，在PC机中，采用基于归一化互信息的配准方法，将图谱向通过CT或核磁扫描得到的个体图像作整体配准，以消除图谱和目标图像间的整体差异；所述图谱是指对已有的医学图象经解剖学专家手工分割后对每个器官区域做好标记的图象，它可以为PC机提供人体解剖结构的参考信息；所述的归一化互信息nMI(S，T*A)用下式表示nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>其中，S为目标图像，T*A代表对图谱A做空间变换T后的图像；H(S)代表目标图像的熵，H(T*A)代表变换后图谱的熵，H(S，T*A)代表目标图像和变换后图谱间的联合熵；nMI是英文的归一化互信息的缩写(normalized mutual information)，nMI(S，T*A)表示S与T*A之间的归一化互信息；空间变换T为包含9个参数的仿射变换T＝(p，q，r，u，v，w，φ，ω，θ)；其中，p，q，r分别为空间变换后x，y，z方向的位移；r，u，v分别为空间变换后x，y，z方向的比例缩放；φ，ω，θ分别为空间变换后绕x，y，z轴的转角；若变换前点坐标为X＝(x，y，z)，则变换后的坐标为X′＝(x′，y′，z′)，x′y′z′1=100p010q001r0001×cosφsinφ00-sinφcosφ0000100001×00-sinω00100sinω0cosω00001]]>×10000cosθsinθ00-sinθcosθ00001×u0000v0000w00001×xyz1]]>至此，定义了目标图像S和经过空间变换T后的图谱A之间的归一化互信息；所述基于归一化互信息的配准方法依次按以下步骤实现步骤(1.1).随机确定初始的空间变换T0＝(p0，q0，r0，u0，v0，w0，φ0，ω0，θ0)中九个参数值，并用由此得到的空间变换作用于图谱A，得到变换后的图谱T0*A；步骤(1.2).计算目标图像S和经过空间变换T后的图谱T*A间的联合归一化联合直方图h(l，k)，h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}其中，S(x)＝l表示在目标图像中像素x对应的灰度值为l，T*A(x)＝k表示在经过空间变换的图谱T*A像素中x对应的灰度值为k；Num为像素x的数目；步骤(1.3).计算目标图像S和经过空间变换T后的图谱T*A间的联合概率分布pS，T*A(l，k)pS,T*A(l,k)=h(l,k)Σl,kh(l,k).]]>步骤(1.4).分别计算图像S的边缘概率分布pS(l)和T*A的边缘概率分布pT*A(k)PS(l)=ΣkPS,T*A(l,k),]]>pT*A(k)=ΣlpS,T*A(l,k).]]>边缘概率分布是概率统计理论中的概念，对于两个随机变量x，y组成的随机向量(x，y)，其分量x的概率分布称为随机向量(x，y)的关于x的边缘分布；这里是把图像S和变换后的图谱T*A分别作为两个随机变量l和k，然后分别计算l和k的边缘概率分布；步骤(1.5).按以下公式分别计算目标图像的熵H(S)，变换后图谱的熵H(T*A)以及目标图像和变换后图谱间的联合熵H(S，T*A)；H(S)=-ΣlpS(l)logpS(l)]]>H(T*A)=-ΣkpT*A(k)logpT*A(k)]]>H(S,T*A)=-ΣkpS,T*A(l,k)logS,T*A(l,k)]]>步骤(1.6).按下式计算目标图像S和经过T变换后的图谱T*A间的归一化互信息nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>步骤(1.7).将归一化互信息nMI(S，T*A)带入Powell优化算法对空间变换T进行反复优化迭代，一直到收敛为止；Powell优化算法是一种传统的确定性优化方法，并且被认为是目前最有效的无需求导数的优方法；该算法的基本思想是，对于n维极值问题，首先沿着n个坐标方向求极小，经过n次之后得到n个共轭方向，然后沿n个共轭方向求极小，经过多次迭代后便可求得极小值；在这里，利用Powell优化算法对空间变换T的值进行优化，以求得到互信息nMI(S，T*A)的极小值；空间变换T含有九个分量，即(p，q，r，u，v，w，φ，ω，θ)，这九个分量分别对应Powell优化算法中的九个维度，即利用Powell优化算法求九维极值问题；当Powell优化算法搜索结束，即算法收敛后，得到的nMI(S，T*A)的极小值所对应的空间变换值即为最优空间变换T*，T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；则T**A为整体配准所要得到的配准图谱；步骤(2).依次按以下步骤对感兴趣的每一个器官进行单独配准，确定每个感兴趣的器官在目标图像S中的大致对应区域；具体分为以下几步步骤(2.1).按下式为每个感兴趣器官确定灰度范围区间RkRk＝[μk-λσk，μk+λσk]，其中，μk，σk为变换后图谱T*A的第k个感兴趣器官在目标图像S中所占据的区域内的灰度均值和标准方差；λ为一个常数，用来调整灰度均匀范围；试验发现，λ的取值范围为1.2～1.3；步骤(2.2).将步骤(1)中得到的最终的空间变换T*作对为每个器官单独进行配准的初始空间变换T′，即令T′＝T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；步骤(2.3).按下式计算变换后的器官区域的相似性测度Mk(T′)Mk(T′)=Nkin(T′)-Nkout(T′).]]>其中，下标k代表对第k个器官进行配准，Nkin(T′)代表采用空间变换T′时器官k在配准图谱中对应区域内包含的器官像素数目，Nkout(T′)代表该对应区域内包含的非器官像素数目；Nkin(T′)=N{x|x∈Ck(T′),S(x)∈Rk}Nkout(T′)=N{x|x∈Ck(T′),S(x)∉Rk}]]>其中，N{}表示计算集合中包含的像素数目，Ck(T′)表示器官k在空间变换T下对应的区域，S(x)表示像素x在CT或MRI图像中的灰度值；步骤(2.4).将变换后的器官区域的相似性测度Mk(T′)带入Powell优化算法对T′进行反复迭代，直到收敛为止；具体过程即利用Powell优化算法对空间变换T′的值进行优化，以求得到相似性测度Mk(T′)的极小值；空间变换T′含有九个分量，即(p′，q′，r′，u′，v′，w′，φ′，ω′，θ′)，这九个分量分别对应Powell优化算法中的九个维度，即利用Powell优化算法求九维极值问题；当Powell优化算法搜索结束，即算法收敛后，得到的Mk(T′)的极小值所对应的空间变换值即为最优空间变换T*′，T′*＝(p′*，q′*，r′*，u′*，v′*，w′*，φ′*，ω′*，θ′*)；则Ck(T′*)为单独配准所要得到的器官区域；步骤(3).模糊连接分割；具体分为以下几步步骤(3.1).设Ck(T′*)为模糊连接对第k个器官进行分割的初始区域，按下式计算该区域中的灰度直方图Hk(l)Hk(l)＝N{x|x∈Ck，S(x)＝l}，l＝0，1……LHk(l)可理解为Ck(T′*)区域中灰度为l的像素x的个数；步骤(3.2).采用曲线拟合的方法对直方图进行平滑和去噪，得到拟合后的直方图 H^k(l)=Σj=1najexp[-(l-bjcj)2],]]>其中，n代表高斯分量的数目，即所要分割的器官的数目；例如我们对肝脏和脾脏两种器官进行分割，则n＝2；aj，bj，cj依次代表第j个高斯分布的峰值、中心和宽度；步骤(3.3).拟合后的直方图曲线可以看作器官在目标图像中的灰度概率密度函数；选择对应最大灰度概率的至少一个像素作为器官的初始种子点；步骤(3.4).考察每一个种子点相邻区域的灰度概率之和；相邻区域定义为与作为中心像素的种子点距离小于一个像素单位的距离的像素集合；所述灰度概率为所研究器官区域内具有该灰度的像素数目与区域内像素总数的比值；具有最大区域灰度概率总和的点将被最终确定为种子点；步骤(3.5).以由步骤(3.2)得到的直方图 所对应的器官区域灰度分布的概率密度函数和最终种子点的位置为初始参数，用下式所述的用模糊连接算法计算整幅图像中每一像素点与器官最终种子点的模糊连接强度cμ(p,q)=maxρ(p,q)∈p(p,q)[minz∈ρ(p,q)H^k(l(z))]]]>其中，p表示最终的种子点，q表示整幅图像中任意一点，cμ(p，q)表示p与q之间的模糊连接强度；ρ(p，q)表示p与q之间的任意一条连通路径，P(p，q)表示p与q之间所有连通路径的集合，z∈ρ(p，q)表示Z是路径ρ(p，q)上的任意一像素点，l(z)为像素Z在目标图像中的灰度值， 表示灰度值l(z)在直方图 中所对应的概率密度值；模糊连接强度表征的该像素属于器官区域的隶属度，将每个像素点的隶属度作为灰度所得到的图像称为模糊连接图像；步骤(3.6).为每个器官的模糊连接图像指定一个最佳阈值；所谓阈值，就是用来将属于目标器官的像素和背景像素区分开的一个数值，规定灰度大于或等于此阈值的像素为属于目标器官的像素，灰度小于此阈值的像素为背景像素；所谓最佳阈值，就是由此阈值得到的器官像素所集合成的区域最符合真实的器官区域；最佳阈值是将模糊连接图像转化为最终的分割结果的重要参量；确定最佳阈值具体步骤如下首先将阈值设定为1，然后逐步降低阈值；用于每一次降低得到的阈值来分割模糊连接图像，得到器官的分割区域；观察分割区域的形状随阈值的改变，如果在某一次阈值降低时，形状剧烈改变，则停止降低阈值，将上一阈值定为最佳阈值；关于分割区域的形状改变的衡量办法和最佳阈值的确定方法，进一步详细解释如下a).分割区域的形状改变的衡量的办法采用面积相对变化值和紧密度相对变化值对来衡量分割区域的形状改变；紧密度(Compactness)的定义为 面积相对变化值和紧密度相对变化值定义如下对于每次阈值降低，分别计算降低前后的器官面积和紧密度；令降低前的器官面积和紧密度为原面积和原紧密度，降低后的器官面积和紧密度为新面积和新紧密度，则面积和紧密度的相对变化定义为 b).用面积相对变化值和紧密度相对变化值确定最优阈值的方法在阈值降低过程中，如果在某一次降低时，发现面积和紧密度的相对变化值大于预先设定的参考值ρ，则认为形状发生了剧烈改变，从而停止阈值搜索，并将降低前的阈值定为最佳阈值；试验表明，参考值ρ取值0.15～0.2能够取得较准确的最佳阈值；步骤(3.7)、按照下式用步骤(3.6)得到的最佳阈值对模糊连接图像进行分割，得到器官的二值分割图像设第k个器官的最佳阈值为Tk，则该器官的二值分割结果为 其中，x代表任一像素，Fk(x)为第k个器官的模糊连接图像在像素x处的灰度值，Bk(x)为分割后的二值图像在像素x处的灰度值；步骤(4).器官形状修正，以解决因为人体腹腔内相互紧贴的两个器官因为灰度相近而造成模糊连接分割可能产生的错误分割结果；医学图像中，经常出现两个在人体腹腔中相互紧贴的器官在图像中的灰度也相似的现象；对于这样的两个器官，医学图像中在其表面相贴处可能没有明显的边界，这会使模糊连接算法认为这两个器官是长在一起的，从而错误地分割成同一个器官；本研究提出以器官形状稳定为依据来检测错误分割，并采用距离变换结合分水岭分割算法的对错误分割结果进行形状修正；具体流程如下步骤(4.1).对于多层医学图像，按以下步骤从上到下，采用逐层分析的方法检测存在错误分割的层面；步骤(4.1.1).设定当前层为算法正在分析的层面，当前层的序号为i下方相邻层为与当前层相邻的下方层面，下方相邻层的序号为i-1；Oki为当前层中的器官分割区域；Oki-1为下方相邻层中器官分割区域；在当前层中，令器官分割区域内的各点到区域边界的距离为负值，对于整幅图像中所有不属于该区域的点，即所有该区域以外的点，令其到区域边界的距离为正值；步骤(4.1.2).对当前层i中的器官分割区域Oki进行距离变换，即对整幅图像中每点按下式计算它们到区域边界的最近距离Dki 其中，y代表了区域边界的像素点中与x距离最近一点；Bik(x)为步骤(3.7)中得到的最终分割结果Bk*(x)在当前层i中的像素值；将距离值Dki(x)作为灰度值赋给每一像素点，就得到了距离图像Dki；步骤(4.1.3).提取下方相邻层中区域Oki-1的边界上各点像素的位置，则在Dki中的相同位置处的各点的Dki(x)值可认为是Oki-1边界上的像素的距离值；得到所有Oki-1的边界像素的距离值，计算这些距离值的标准方差，以确定器官在第i层和第i-1层分割结果间的相似程度；如果标准方差大于预先设定的参考值，则认为检测到了分割错误，下一层即为存在错误分割的层面；根据经验，这个参考值一般设定在5个像素距离；按照上述方法，从上到下检测所有层面，得到所有存在错误分割的层面；步骤(4.2).按下式对存在错误分割的各层面的二值图像进行距离变换；设定对于器官区域内所有像素点，令其距离值为负值；对于整幅图像中区域外所有不属于该区域的区域以外的点，令其距离值为负的极小值，即负无穷-∞；则距离图像Dki′(x) 其中y代表了与x距离最近的边界点；Bki即存在错误分割的层面的二值图像；步骤(4.3).按以下方法对步骤(4.2)得到的距离图像Dki′应用分水岭算法；步骤(4.3.1).对整幅图像设定一个灰度阈值T；令T的初值为图像Dki′中除了负无穷以外的最小值，即T0=min({Dki(x)|Dki(x)≠-∞}),]]>其中，x代表图像中任一个像素点的坐标；步骤(4.3.2).以d为步长逐渐增大阈值T，即Tn＝T0+n·d，其中，Tn表示第n次增加步长时阈值T的取值； n取值范围1≤n≤K，其中K=(max(Dki(x))-T0)/d]]>，max(Dki(x))表示图像Dki′的最大值；d取值固定，根据经验取d＝3；步骤(4.3.3).利用第n次增加步长得到的阈值Tn按下式将图像Dki′分割成二值图像Wn 其中，x代表图像中任一个像素点的坐标；Wn(x)＝1表示淹没在水面以下的部分，Wn(x)＝0表示露在水面以上的陆地部分；对于前一次(即第n-1次)增加步长得到阈值Tn-1，按下式将图像Dki′分割成二值图像Wn-1 其中，x代表图像中任一个像素点的坐标；步骤(4.3.4).在Tn的增长过程中，每当Tn大于一个局部极大值，就说明图像Wn-1内有两个连通区域在Wn中合并成了一个连通域，则将Wn-1和Wn两幅图像进行异或操作，得到分水岭图像Sn 其中，x代表图像中任一个像素点的坐标，表示异或操作，Sn(x)＝1表示像素x为分水岭所在位置；步骤(4.3.5).对n取值范围1≤n≤K中的每一个n值，均会得到一个相应的分水岭图像Sn；设所有的n值所对应的分水岭图像为S1，S2，…，Smax(Dki(x))，则分水岭变换的结果图像S为所有分水岭图像的集合S=S1·S2·...·Smax(Dki(x));]]>图像S中所有像素值为1的像素就表示图像Dki′中的分水岭所在的位置；由这些分水岭即把图像Dki′分割成若干个子区域；步骤(4.4).对分水岭算法得到的图像Dki′中的每个子区域，分别计算其与下方相邻层面器官分割区域Oki-1的面积重合比率 把重合比率低于50％的区域视为错误分割区域，重新标记为背景区域，从而将其在距离图像Dki′中除去。
全文摘要
本发明属于医学图像中腹部器官分割技术领域。其特征在于，该方法依次含有以下步骤采用基于归一化互信息的配准方法将标准图谱向通过CT或核磁扫描得到的个体图像作整体配准；用相似性测度对每个感兴趣器官进行单独配准；模糊连接分割和器官形状修正；其中，模糊连接分割又分为先算出曲线拟和后的灰度直方图曲线，并视其为目标图像中灰度概率密度的函数；再把具有最大区域灰度概率密度的点作为种子点，将直方图拟和曲线以及种子作为初始参数，用模糊连接算法得到模糊连接图像，再用设定的阈值把目标像素和背景像素分开。本发明具有可同时进行目标图像中腹部多器官自动分割，实施简单的优点。
文档编号A61B5/055GK101013503SQ20071006306
公开日2007年8月8日 申请日期2007年1月26日 优先权日2007年1月26日
发明者白净, 周永新, 王洪凯, 刘加成, 张永红, 张菊鹏 申请人:清华大学