本发明提供一种先验形状约束的图像水平集分割方法，涉及计算机视觉及图像处理技术领域。该方法读取待分割图像，并初始化图像的水平基函数和形状先验水平集函数；建立图像水平集分割的能量泛函，设定水平集方法中的参数；然后分别更新第j次迭代的聚簇中心集合、图像水平集函数、形状约束水平集函数和基函数权重列向量，并计算第j次迭代的能量泛函的函数值直至达到最大迭代次数J或者能量泛函相邻两次的迭代结果差值小于设定的收敛阈值，得到待分割图像I的分割结果，并根据更新后的基函数权重列向量和基函数列向量得到待分割图像的偏移场估计；本发明方法克服了图像中灰度不均对分割精度造成的影响，在先验形状的约束下，增强了分割方法的鲁棒性。
1.一种先验形状约束的图像水平集分割方法，其特征在于：包括以下步骤：步骤1：读取待分割图像I作为目标图像；
步骤2：选定正交基函数集中的前M个低阶函数，用于在目标图像的定义域Ω内拟合图像灰度的不一致性；
步骤3：初始化划分目标图像区域的水平集函数：设定初始的零水平集分割曲线，将目0
标图像划分成2个区域，并根据设定结果生成高一维水平集函数φ (x)，其中，x为定义域Ω内任意一像素点；
步骤4：初始化标记待划分目标图像区域形状特征的形状约束水平集函数：构造符合感兴趣区域形状特征的闭合曲线，将其作为形状约束水平集函数的零水平集曲线，并根据设定的零水平集曲线生成高一维的形状约束水平集函数步骤5：建立目标图像水平集分割的能量泛函 其中，c＝(c1,c2)为聚簇中心T集合，c1，c2是零水平集分割曲线划分成的2个子区域内的灰度均值，w＝(w1,w2,...,wM) 为由基函数权重值w1,w2,...,wM构成的列向量；并设定图像水平集分割方法中的参数；
步骤6：给定选定基函数的初始权重值，即 并设定最大迭代次数J和目标图像水平集分割能量泛函 的收敛阈值ε；
步骤7：分别更新第j次迭代的聚簇中心集合c、图像水平集函数φ 、形状约束水平集函j数 基函数权重列向量w ，并依据水平集函数值判别当前像素点x归属于划分区域1或2的j j j j可能性度量，即隶属度函数M1(φ (x))＝H(φ (x)),M2(φ (x))＝1‑H(φ (x))，并计算第j次迭代的能量泛函 函数值；
步骤8：判断是否达到迭代终止条件：若当前迭代次数j达到最大迭代次数J或者能量泛函 相邻两次的迭代结果差值小于设定的收敛阈值ε，则执行步骤9，否则，迭代次数j加1，返回步骤7；
步骤9：根据当前更新后的图像水平集函数φ (x)构造图像的隶属度函数M1(φ (x)),M2j j(φ (x))，即待分割图像I的分割结果，并根据更新后的基函数权重列向量w 和基函数列向量G得到待分割图像I的偏移场估计b；
步骤1所述图像I定义域为二维实数空间子集Ω，给定该定义域内任意一像素点x，受噪声与灰度不一致问题的影响，目标图像在该点的灰度值表示为：I(x)＝b(x)R(x)+n(x)其中，R(x)为像素点x的理想灰度值，其仅与成像协议以及待成像物体质地有关，b(x)为灰度不一致性问题在像素点x的偏移场，n(x)为噪声在x点的取值；
所述灰度不一致性问题在像素点x的偏移场b(x)如下公式所示：T
其中，w＝(w1,w2,...,wM) 为由基函数权重值w1,w2,...,wM构成的列向量，G(x)＝(g1T(x),g2(x),...,gM(x)) 为两两正交的基函数g1,g2,...,gM在像素点x处的取值构成的列向量；
步骤5所述能量泛函如下公式所示：E＝αED+βES+γEO
其中，ED为目标图像数据项，ES为感兴趣区域形状约束项，EO为正则项，x与y均为定义域为图像域的像素变量，H为Heaviside函数；
所述设定的图像水平集分割方法中的参数为：加权系数α、加权系数β、加权系数γ、水平集函数正则项系数μ、水平集分割曲线平滑项系数v及平移因子a；
步骤7所述更新聚簇中心集合c、形状约束水平集 图像水平集函数φ 、基函数权重列j向量w的具体方法为：
所述更新聚簇中心集合c的公式如下所示：其中，下标n表示划分区域编号，其取值为1或2；
所述更新形状约束水平集函数 由初始的形状水平集函数 通过公式变换得到；
其中， a 、r和θ分别为第j次迭代的平移因子、放缩因子和旋转因子；
所述更新图像水平集函数φ的公式如下所示：其中，Δt为时间步长， 为迭代更新的水平集函数的时间偏导数；所述时间偏导数公式如下所示：其中，δ是狄拉克函数，其为heaviside函数的导数；
其中，G(y)＝(g1(y),g2(y),...,gM(y)) 。
2.根据权利要求1所述的一种先验形状约束的图像水平集分割方法，其特征在于：所述设定的初始零水平集分割曲线为指定形状的闭合曲线或随机生成的闭合曲线；所述水平集0
函数φ (x)为根据设定的初始零水平集分割曲线生成的高一维函数阶梯函数或符号距离函数；所述形状约束水平集函数 为根据设定的初始零水平集曲线生成的高一维函数阶梯函数或符号距离函数。
3.根据权利要求2所述的一种先验形状约束的图像水平集分割方法，其特征在于：步骤j j j
9所述根据当前更新后的图像水平集函数φ (x)构造图像的隶属度函数M1(φ (x)),M2(φ(x))由水平集函数构造得出，且若待分割图像I中的像素点x属于待划分的2个区域中的第nj j个子区域Ωn中，则隶属度函数Mn(φ(x))＝1，否则Mn(φ(x))＝0。
一种先验形状约束的图像水平集分割方法
技术领域
[0001] 本发明涉及计算机视觉及图像处理技术领域，尤其涉及一种先验形状约束的图像水平集分割方法。
背景技术
[0002] 图像分割在计算机视觉和图像处理领域起重要作用，目的是将图像域中感兴趣的目标区域划分出来，为后续处理分析特别是图像理解提供技术支撑。水平集分割方法用平滑封闭曲线标记感兴趣目标区域边界，通过不断迭代演化实现分割图像的目的。该类方法易进行维度扩展，能够提供亚像素级的分割精度，被广泛应用于图像分割领域。但是，传统水平集方法只考虑了感兴趣目标区域与背景区域间的灰度差异，并未充分利用感兴趣目标区域的更高维度特征。在诸多应用场景中，如图1所示，图像中感兴趣目标区域具有特定形状特点。为使水平集方法能够准确从图像中提取具有特定形状特征的感兴趣目标区域，需要一种先验形状约束的图像水平集分割方法。
发明内容
[0003] 本发明要解决的技术问题是针对上述现有技术的不足，提供一种先验形状约束的图像水平集分割方法，实现对感兴趣区域具有特定形状特征的图像进行水平集分割的目的。
[0004] 为解决上述技术问题，本发明所采取的技术方案是：一种先验形状约束的图像水平集分割方法，包括以下步骤：
[0005] 步骤1：读取待分割图像I作为目标图像；图像I定义域为二维实数空间子集Ω，给定该定义域内任意一像素点x，受噪声与灰度不一致问题的影响，目标图像在该点的灰度值表示为：
[0006] I(x)＝b(x)R(x)+n(x)
[0007] 其中，R(x)为像素点x的理想灰度值，其仅与成像协议以及待成像物体质地有关，b(x)为灰度不一致性问题在像素点x的偏移场，n(x)为噪声在x点的取值；
[0008] 步骤2：选定正交基函数集中的前M个低阶函数，用于在目标图像的定义域Ω内拟合图像灰度的不一致性；具体地，所述灰度不一致性问题在像素点x的偏移场b(x)如下公式所示：
[0009] b(x)＝wTG(x)
[0010] 其中，w＝(w1，w2，...，wM)T为由基函数权重值w1，w2，...，wM构成的列向量，G(x)＝T(g1(x)，g2(x)，...，gM(x)) 为两两正交的基函数g1，g2，...，gM在像素点x处的取值构成的列向量；
[0011] 步骤3：初始化划分目标图像区域的水平集函数：设定初始的零水平集分割曲线，0
将目标图像划分成2个区域，并根据设定结果生成高一维水平集函数φ(x)；
[0012] 所述设定的初始零水平集分割曲线为指定形状的闭合曲线或随机生成的闭合曲0
线；所述水平集函数φ (x)为根据设定的初始零水平集分割曲线生成的高一维函数阶梯函数或符号距离函数；
[0013] 步骤4：初始化标记待划分目标图像区域形状特征的形状约束水平集函数：构造符合感兴趣区域形状特征的闭合曲线，将其作为形状约束水平集函数的零水平集曲线，并根据设定的零水平集曲线生成高一维的形状约束水平集函数
[0014] 所述形状约束水平集函数 为根据设定的初始零水平集曲线生成的高一维函数阶梯函数或符号距离函数；
[0015] 步骤5：建立目标图像水平集分割的能量泛函 其中，c＝(c1，c2)为聚簇中心集合，c1，c2是零水平集分割曲线划分成的2个子区域内的灰度均值；并设定图像水平集分割方法中的参数：加权系数α、加权系数β、加权系数γ、水平集函数正则项系数μ、水平集分割曲线平滑项系数v及平移因子a；
[0016] 所述能量泛函如下公式所示：
[0017] E＝αED+βES+γEO
[0018] 其中，
[0019]
[0020]
[0021]
[0022] 其中，ED为目标图像数据项，ES为感兴趣区域形状约束项，EO为正则项，x与y均为定义域为图像域的像素变量，H为Heaviside函数；
[0023] 步骤6：给定选定基函数的初始权重值，即 并设定最大迭代次数J和目标图像水平集分割能量泛函 的收敛阈值ε；
[0024] 步骤7：分别更新第j次迭代的聚簇中心集合cj、图像水平集函数φj、形状约束水平j集函数 基函数权重列向量w ，并依据水平集函数值判别当前像素点x归属于划分区域1j j j j
或2的可能性度量，即隶属度函数M1(φ (x))＝H(φ (x))，M2(φ (x))＝1‑H(φ (x))，并计算第j次迭代的能量泛函 函数值；
[0025] 所述更新聚簇中心集合cj、形状约束水平集 图像水平集函数φj、基函数权重列j向量w的具体方法为：
[0026] 所述更新聚簇中心集合cj的公式如下所示：
[0027]
[0028] 其中，下标n表示划分区域编号，其取值为1或2；
[0029] 所述更新形状约束水平集函数 由初始的形状水平集函数 通过公式变换得到；
j j j
[0030] 其中， a 、r 和θ分别为第j次迭代的平移因子、放缩因子和旋转因子；
[0031] 平移因子aj的更新公式如下所示：
[0032]
[0033] 放缩因子rj的更新公式如下所示：
[0034]
[0035] 旋转因子θj的更新公式如下所示：
[0036]
[0037] 所述更新图像水平集函数φj的公式如下所示：
[0038]
[0039] 其中，Δt为时间步长， 为迭代更新的水平集函数的时间偏导数；所述时间偏导数公式如下所示：
[0040]
[0041] 其中，δ是狄拉克函数，其为heaviside函数的导数；
[0042] 所述更新基函数权重列向量wj的公式如下所示：
[0043]
[0044] 其中，G(y)＝(g1(y)，g2(y)，...，gM(y))T；
[0045] 步骤8：判断是否达到迭代终止条件：若当前迭代次数j达到最大迭代次数J或者能量泛函 相邻两次的迭代结果差值小于设定的收敛阈值ε，则执行步骤9，否则，迭代次数j加1，返回步骤7；
[0046] 步骤9：根据当前更新后的图像水平集函数φj(x)构造图像的隶属度函数M1(φjj j(x))，M2(φ (x))，即待分割图像I的分割结果，并根据更新后的基函数权重列向量w 和基函数列向量G得到待分割图像I的偏移场估计b；
[0047] 所述根据当前更新后的图像水平集函数φj(x)构造图像的隶属度函数M1(φjj(x))，M2(φ (x))由水平集函数构造得出，且若待分割图像I中的像素点x属于待划分的2个j j
区域中的第n个子区域Ωn中，则隶属度函数Mn(φ(x))＝1，否则Mn(φ(x))＝0。
[0048] 采用上述技术方案所产生的有益效果在于：本发明提供的一种先验形状约束的图像水平集分割方法，克服了图像中的灰度不均对分割精度造成的影响，在先验形状的约束下，增强了分割方法的鲁棒性。应用该方法能够得到精确的图像分割结果，同时估计和校正了图像的偏移场。
附图说明
[0049] 图1为本发明背景技术提供的左心室短轴核磁图像；
[0050] 图2为本发明实施例1提供的一种先验形状约束的图像水平集分割方法的流程图；
[0051] 图3为本发明实施例1提供的设定的初始零水平集分割曲线和形状水平集函数零水平集曲线，其中，(a)为设定的初始零水平集分割曲线，(b)为设定的初始形状水平集函数零水平集曲线；
[0052] 图4为本发明实施例1提供的水平集函数分割结果及金标准分割结果，其中，(a)为迭代10次后的水平集函数分割结果，(b)为迭代30次后的水平集函数分割结果，(c)为迭代50次后的水平集函数分割结果，(d)为实施例1中的金标准分割结果；
[0053] 图5为本发明实施例1得到的待分割图像I的偏移场估计及灰度不一致性得到校正的图像，其中，(a)为得到的待分割图像I的偏移场估计，(b)为灰度不一致校正后的图像；
[0054] 图6为本发明实施例2提供的冠状主动脉图像和冠状主动脉分割结果，其中，(a)为待分割的冠状主动脉图像和初始零水平集分割曲线，(b)为待分割图像的偏移估计图像，(c)为偏移场校正后的图像，(d)为(a)中冠状动脉的分割结果。
具体实施方式
[0055] 下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
[0056] 实施例1：
[0057] 一种先验形状约束的图像水平集分割方法，如图2所示，包括以下步骤：
[0058] 步骤1：读取待分割图像I作为目标图像；图像I定义域为二维实数空间子集Ω，给定该定义域内任意一像素点x，受噪声与灰度不一致问题的影响，目标图像在该点的灰度值表示为：
[0059] I(x)＝b(x)R(x)+n(x)
[0060] 其中，R(x)为像素点x的理想灰度值，其仅与成像协议以及待成像物体质地有关，b(x)为灰度不一致性问题在像素点x的偏移场，n(x)为噪声在x点的取值；
[0061] 步骤2：选定正交基函数集中的前M个低阶函数，用于在目标图像的定义域Ω内拟合图像灰度的不一致性；具体地，所述灰度不一致性问题在像素点x的偏移场b(x)如下公式所示：
[0062] b(x)＝wTG(x)
[0063] 其中，w＝(w1，w2，...，wM)T为由基函数权重值w1，w2，...，wM构成的列向量，G(x)＝T(g1(x)，g2(x)，...，gM(x)) 为两两正交的基函数g1，g2，...，gM在像素点x处的取值构成的列向量；
[0064] 步骤3：初始化划分目标图像区域的水平集函数：设定初始的零水平集分割曲线，0
将目标图像划分成2个区域，并根据设定结果生成高一维水平集函数φ(x)；
[0065] 所述设定的初始零水平集分割曲线为指定形状的闭合曲线或随机生成的闭合曲0
线；所述水平集函数φ (x)为根据设定的初始零水平集分割曲线生成的高一维函数阶梯函数或符号距离函数；
[0066] 步骤4：初始化标记待划分目标图像区域形状特征的形状约束水平集函数：构造符合感兴趣区域形状特征的闭合曲线，将其作为形状约束水平集函数的零水平集曲线，并根据设定的零水平集曲线生成高一维的形状约束水平集函数
[0067] 所述形状约束水平集函数 为根据设定的初始零水平集曲线生成的高一维函数阶梯函数或符号距离函数；
[0068] 本实施例中，目的是分割出短轴核磁图像中左心室的区域，因此待划分的区域数为2，设定的初始零水平集分割曲线为闭合圆。同样地，将形状水平集函数 的零水平集初始为闭合圆形曲线，如图3所示。
[0069] 步骤5：建立目标图像水平集分割的能量泛函 其中，c＝(c1，c2)为聚簇中心集合，c1，c2是零水平集分割曲线划分成的2个子区域内的灰度均值；并设定图像水平集分割方法中的参数：加权系数α、加权系数β、加权系数γ、水平集函数正则项系数μ、水平集分割曲线平滑项系数v及平移因子a；
[0070] 所述能量泛函如下公式所示：
[0071] E＝αED+βES+γEO
[0072] 其中，
[0073]
[0074]
[0075]
[0076] 其中，ED为目标图像数据项，ES为感兴趣区域形状约束项，EO为正则项，x与y均为定义域为图像域的像素变量，H为Heaviside函数；
[0077] 本实施例中，图像数据项ED用于将零水平集分割演化曲线吸引到待分割对象边界位置，形状约束项ES用于将零水平集分割曲线φ(x)保持形状水平集函数 相同的形状。
[0078] 本实施例中，图像数据项系数越大，零水平集分割曲线越接近待分割对象在图像中的边界，因此根据待分割图像设定图像水平集分割的控制参数：加权系数α＝0.1，加权系数β＝0.9，加权系数γ＝1，水平集函数正则项系数μ＝1.0，水平集分割曲线平滑项系数v＝0.005*255*255，平移因子a＝0。
[0079] 步骤6：给定选定基函数的初始权重值，即 并设定最大迭代次数J和目标图像水平集分割能量泛函 的收敛阈值ε；
[0080] 本实施例中，设定最大迭代次数J为50，图像水平集分割的能量泛函 的收敛阈值为0.01。
[0081] 步骤7：分别更新第j次迭代的聚簇中心集合cj、图像水平集函数φj、形状约束水平j集函数 基函数权重列向量w ，并依据水平集函数值判别当前像素点x归属于划分区域1j j j j
或2的可能性度量，即隶属度函数M1(φ (x))＝H(φ (x))，M2(φ (x))＝1‑H(φ (x))，并计j
算第j次迭代的能量泛函 函数值；所述更新聚簇中心集合c 、形状约束水平集j j
图像水平集函数φ、基函数权重列向量w的具体方法为：
[0082] 本实施方式中，当零水平集分割曲线正好位于心脏核磁图像中左心室边界时，公式(3)给出的能量泛函 取得最小值，因此，为使得 取得最小值，保持j
φ， 与w不变，解 即更新聚簇中心集合c的公式如式(7)所示：
[0083]
[0084] 其中，下标n表示划分区域编号，其取值为1或2；
[0085] 所述更新形状约束水平集函数 由初始的形状水平集函数 通过公式变换得到；
[0086] 其中， aj、rj和θj分别为第j次迭代的平移因子、放缩因子和旋转因子，r初始为1、θ初始为0，且r和θ始终保持初始值不变；
j
[0087] 本实施例中，保持φ、w与c不变，解 等价于 平移因子a的更新公式如下所示：
[0088]
[0089] 放缩因子rj的更新公式如下所示：
[0090]
[0091] 旋转因子θj的更新公式如下所示：
[0092]
[0093] 本实施方式中，保持w、c与 不变，采用最大梯度降方法求得更新后的图像水平集j函数φ的公式如下所示：
[0094]
[0095] 其中，Δt为时间步长， 为迭代更新的水平集函数的时间偏导数；所述时间偏导数公式如下所示：
[0096]
[0097] 其中，δ是狄拉克函数，其为heaviside函数的导数；
[0098] 本实施例中，保持φ、 与c不变，解 得到更新后的基函数权重列向量jw的公式如下所示：
[0099]
[0100] 其中，G(y)＝(g1(y)，g2(y)，...，gM(y))T；
[0101] 步骤8：判断是否达到迭代终止条件：若当前迭代次数j达到最大迭代次数J或者能量泛函 相邻两次的迭代结果差值小于设定的收敛阈值ε，则执行步骤9，否则，迭代次数j加1，返回步骤7；
[0102] 本实施例中，迭代10次、30次和50次后的水平集函数分割结果如图4所示。
[0103] 步骤9：根据当前更新后的图像水平集函数φj(x)构造图像的隶属度函数M1(φjj j(x))，M2(φ (x))，即待分割图像I的分割结果，并根据更新后的基函数权重列向量w 和基函数列向量G得到待分割图像I的偏移场估计b；
[0104] 所述根据当前更新后的图像水平集函数φj(x)构造图像的隶属度函数M1(φjj(x))，M2(φ (x))由水平集函数构造得出，且若待分割图像I中的像素点x属于待划分的2个j j
区域中的第n个子区域Ωn中，则隶属度函数Mn(φ(x))＝1，否则Mn(φ(x))＝0。
[0105] 本实施例中，M1(φj(x))，M2(φj(x))由水平集函数的Heaviside变换的逻辑运算得出，且待分割心脏核磁图像I中的像素点x属于待划分的2个区域中的第n个子区域Ωn中，j j则隶属度函数Mn(φ(x))＝1，否则Mn(φ(x))＝1。
[0106] 本实施例中，得到待分割图像I的偏移场估计、灰度不一致性得到校正的图像如图5所示。
[0107] 实施例2：
[0108] 本实施对存在噪声、偏移场的冠状动脉图像采用本发明方法进行分割，利用冠状主动脉的形状特点作为先验形状约束，将图像分割为两个区域，如图6(a)所示，本发明方法得到的待分割图像的偏移场估计如图6(b)所示，偏移场校正后的图像如图6(c)所示，得到的待分割图像的最终分割结果如图6(d)所示。
[0109] 最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。
