[发明专利]一种医学图像中腹部器官分割方法

摘要

本发明属于医学图像中腹部器官分割技术领域。其特征在于，该方法依次含有以下步骤：采用基于归一化互信息的配准方法将标准图谱向通过CT或核磁扫描得到的个体图像作整体配准；用相似性测度对每个感兴趣器官进行单独配准；模糊连接分割和器官形状修正；其中，模糊连接分割又分为：先算出曲线拟和后的灰度直方图曲线，并视其为目标图像中灰度概率密度的函数；再把具有最大区域灰度概率密度的点作为种子点，将直方图拟和曲线以及种子作为初始参数，用模糊连接算法得到模糊连接图像，再用设定的阈值把目标像素和背景像素分开。本发明具有可同时进行目标图像中腹部多器官自动分割，实施简单的优点。

主权项

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).随机确定初始的空间变换T⁰＝(p⁰，q⁰，r⁰，u⁰，v⁰，w⁰，φ⁰，ω⁰，θ⁰)中九个参数值，并用由此得到的空间变换作用于图谱A，得到变换后的图谱T⁰*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间的联合概率分布p_S，T*A(l，k)：pS,T*A(l,k)=h(l,k)Σl,kh(l,k).]]>步骤(1.4).分别计算图像S的边缘概率分布p_S(l)和T*A的边缘概率分布p_T*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).按下式为每个感兴趣器官确定灰度范围区间R_kR_k＝[μ_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).按下式计算变换后的器官区域的相似性测度M_k(T′)：Mk(T′)=Nkin(T′)-Nkout(T′).]]>其中，下标k代表对第k个器官进行配准，N_kⁱⁿ(T′)代表采用空间变换T′时器官k在配准图谱中对应区域内包含的器官像素数目，N_k^out(T′)代表该对应区域内包含的非器官像素数目；Nkin(T′)=N{x|x∈Ck(T′),S(x)∈Rk}Nkout(T′)=N{x|x∈Ck(T′),S(x)∉Rk}]]>其中，N{}表示计算集合中包含的像素数目，C_k(T′)表示器官k在空间变换T下对应的区域，S(x)表示像素x在CT或MRI图像中的灰度值；步骤(2.4).将变换后的器官区域的相似性测度M_k(T′)带入Powell优化算法对T′进行反复迭代，直到收敛为止；具体过程即利用Powell优化算法对空间变换T′的值进行优化，以求得到相似性测度M_k(T′)的极小值；空间变换T′含有九个分量，即(p′，q′，r′，u′，v′，w′，φ′，ω′，θ′)，这九个分量分别对应Powell优化算法中的九个维度，即利用Powell优化算法求九维极值问题；当Powell优化算法搜索结束，即算法收敛后，得到的M_k(T′)的极小值所对应的空间变换值即为最优空间变换T^*′，T′^*＝(p′^*，q′^*，r′^*，u′^*，v′^*，w′^*，φ′^*，ω′^*，θ′^*)；则C_k(T′^*)为单独配准所要得到的器官区域；步骤(3).模糊连接分割；具体分为以下几步：步骤(3.1).设C_k(T′^*)为模糊连接对第k个器官进行分割的初始区域，按下式计算该区域中的灰度直方图H_k(l)：H_k(l)＝N{x|x∈C_k，S(x)＝l}，l＝0，1……LH_k(l)可理解为C_k(T′^*)区域中灰度为l的像素x的个数；步骤(3.2).采用曲线拟合的方法对直方图进行平滑和去噪，得到拟合后的直方图H^k(l)=Σj=1najexp[-(l-bjcj)2],]]>其中，n代表高斯分量的数目，即所要分割的器官的数目；例如我们对肝脏和脾脏两种器官进行分割，则n＝2；a_j，b_j，c_j依次代表第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个器官的最佳阈值为T_k，则该器官的二值分割结果为其中，x代表任一像素，F_k(x)为第k个器官的模糊连接图像在像素x处的灰度值，B_k(x)为分割后的二值图像在像素x处的灰度值；步骤(4).器官形状修正，以解决因为人体腹腔内相互紧贴的两个器官因为灰度相近而造成模糊连接分割可能产生的错误分割结果；医学图像中，经常出现两个在人体腹腔中相互紧贴的器官在图像中的灰度也相似的现象；对于这样的两个器官，医学图像中在其表面相贴处可能没有明显的边界，这会使模糊连接算法认为这两个器官是长在一起的，从而错误地分割成同一个器官；本研究提出以器官形状稳定为依据来检测错误分割，并采用距离变换结合分水岭分割算法的对错误分割结果进行形状修正；具体流程如下：步骤(4.1).对于多层医学图像，按以下步骤从上到下，采用逐层分析的方法检测存在错误分割的层面；步骤(4.1.1).设定：当前层为算法正在分析的层面，当前层的序号为i下方相邻层为与当前层相邻的下方层面，下方相邻层的序号为i-1；O_kⁱ为当前层中的器官分割区域；O_k^i-1为下方相邻层中器官分割区域；在当前层中，令器官分割区域内的各点到区域边界的距离为负值，对于整幅图像中所有不属于该区域的点，即所有该区域以外的点，令其到区域边界的距离为正值；步骤(4.1.2).对当前层i中的器官分割区域O_kⁱ进行距离变换，即对整幅图像中每点按下式计算它们到区域边界的最近距离D_kⁱ：其中，y代表了区域边界的像素点中与x距离最近一点；Bⁱ_k(x)为步骤(3.7)中得到的最终分割结果B_k^*(x)在当前层i中的像素值；将距离值D_kⁱ(x)作为灰度值赋给每一像素点，就得到了距离图像D_kⁱ；步骤(4.1.3).提取下方相邻层中区域O_k^i-1的边界上各点像素的位置，则在D_kⁱ中的相同位置处的各点的D_kⁱ(x)值可认为是O_k^i-1边界上的像素的距离值；得到所有O_k^i-1的边界像素的距离值，计算这些距离值的标准方差，以确定器官在第i层和第i-1层分割结果间的相似程度；如果标准方差大于预先设定的参考值，则认为检测到了分割错误，下一层即为存在错误分割的层面；根据经验，这个参考值一般设定在5个像素距离；按照上述方法，从上到下检测所有层面，得到所有存在错误分割的层面；步骤(4.2).按下式对存在错误分割的各层面的二值图像进行距离变换；设定：对于器官区域内所有像素点，令其距离值为负值；对于整幅图像中区域外所有不属于该区域的区域以外的点，令其距离值为负的极小值，即负无穷-∞；则距离图像D_kⁱ′(x)其中y代表了与x距离最近的边界点；B_kⁱ即存在错误分割的层面的二值图像；步骤(4.3).按以下方法对步骤(4.2)得到的距离图像D_kⁱ′应用分水岭算法；步骤(4.3.1).对整幅图像设定一个灰度阈值T；令T的初值为图像D_kⁱ′中除了负无穷以外的最小值，即：T0=min({Dki(x)|Dki(x)≠-∞}),]]>其中，x代表图像中任一个像素点的坐标；步骤(4.3.2).以d为步长逐渐增大阈值T，即T_n＝T₀+n·d，其中，T_n表示第n次增加步长时阈值T的取值；n取值范围1≤n≤K，其中K=(max(Dki(x))-T0)/d]]>，max(D_kⁱ(x))表示图像D_kⁱ′的最大值；d取值固定，根据经验取d＝3；步骤(4.3.3).利用第n次增加步长得到的阈值T_n按下式将图像D_kⁱ′分割成二值图像W_n：其中，x代表图像中任一个像素点的坐标；W_n(x)＝1表示淹没在水面以下的部分，W_n(x)＝0表示露在水面以上的陆地部分；对于前一次(即第n-1次)增加步长得到阈值T_n-1，按下式将图像D_kⁱ′分割成二值图像W_n-1：其中，x代表图像中任一个像素点的坐标；步骤(4.3.4).在T_n的增长过程中，每当T_n大于一个局部极大值，就说明图像W_n-1内有两个连通区域在W_n中合并成了一个连通域，则将W_n-1和W_n两幅图像进行异或操作，得到分水岭图像S_n：其中，x代表图像中任一个像素点的坐标，表示异或操作，S_n(x)＝1表示像素x为分水岭所在位置；步骤(4.3.5).对n取值范围1≤n≤K中的每一个n值，均会得到一个相应的分水岭图像S_n；设所有的n值所对应的分水岭图像为S₁，S₂，…，S_max(Dki(x))，则分水岭变换的结果图像S为所有分水岭图像的集合：S=S1·S2·...·Smax(Dki(x));]]>图像S中所有像素值为1的像素就表示图像D_kⁱ′中的分水岭所在的位置；由这些分水岭即把图像D_kⁱ′分割成若干个子区域；步骤(4.4).对分水岭算法得到的图像D_kⁱ′中的每个子区域，分别计算其与下方相邻层面器官分割区域O_k^i-1的面积重合比率：把重合比率低于50％的区域视为错误分割区域，重新标记为背景区域，从而将其在距离图像D_kⁱ′中除去。