本发明属于图像处理技术领域，尤其涉及一种独立成分分析多形状先验水平集方法、图像分割系统。
背景技术：
在过去的半个世纪中，图像分割一直受到人们的关注并保持持久的研究热度，至今已提出上千种基于不同理论的分割方法。其中基于主动轮廓模型的方法由于在分割过程中采用先验知识进行指导，并且为分割提供了统一的框架，是目前研究很热门的一类分割方法。从根本上讲，主动轮廓模型方法可以分为参数型和几何型两类，参数型是指直接以参数形式来表达形变的曲线或曲面；几何型是将低维曲线或曲面嵌入到更高一维曲面中，作为其零水平集，并通过几何测度(如单位法向、曲度)来描述，具有可以方便地处理拓扑变换等优点，应用广泛。在图像存在噪声，背景较复杂的情况下，只利用图像数据，难以准确分割目标，因此需要引入高层先验知识。目前引入形状先验的方法有广泛的研究，可以分为两大类：基于单先验的方法和基于多先验的方法。基于单先验的方法中只引入一个形状先验，基本思路是定义当前水平集函数与经过仿射变换的形状先验之间的距离，将形状先验引入能量泛函从而直接作用于分割过程。为了排除其他目标的影响，Daniel Cremers等人通过加入动态标签函数以削弱图像中其他目标对曲线演化的影响，实现对多目标的有选择分割。此外，Junyan Wang等人通过首先建立先验形状内部特征点和形状边缘之间的联系，然后将先验形状和待分割图像进行特征点匹配，估计待分割目标的形状边缘，再进一步优化在特征匹配的约束下的主动轮廓模型能量函数，实现了一种无需初始化的自动的单先验分割方法。该类方法的缺点是单个形状先验不够灵活，不能反映一类形状的共性，而待分割目标和形状先验一般不能只用仿射变换(平移、旋转和缩放)进行匹配，所以当待分割目标和形状先验的差异比较大时，不能得到理想的分割结果。基于多先验的方法中，Michael E.Leventon等人使用PCA提取多个形状先验的局部变化规律，然后基于最大后验概率设计形状约束的能量项，结合图像边缘信息实现分割，该方法假设形状先验符合高斯分布，但实际中往往不能满足该要求。Samuel Dambreville等人使用Kernel-PCA代替上述Leventon文章中使用的PCA，以便于在子空间中更准确地捕获形状先验的统计规律。Daniel Cremers等人使用Parazen窗建立统计模型并构造形状驱动的能量项，将形状信息和图像区域信息结合起来进行图像分割，该方法更加符合实际情况。但是，目前上述多先验方法存在的问题是：由于一般情况下先验形状数目远远小于特征维数，因此先验形状稀疏地分布在高维特征空间，难以估计准确的分布。因此当所估计的形状先验的分布不准确时，形状约束能量项不能有效地约束曲线演化，从而导致分割结果不准确。
综上所述，现有技术存在的问题是：目前引入形状先验的方法中，对于单先验方法，存在单个形状先验不能反映一类形状的共性，当待分割目标和形状先验差异比较大时，不能得到较好的效果的问题；对于多形状先验的方法，存在先验形状一般稀疏地分布在高维特征空间，难以准确地估计形状先验的分布，不准确的分布估计会导致不准确的分割结果的问题。
技术实现要素：
针对现有技术存在的问题，本发明提供了一种独立成分分析多形状先验水平集方法、图像分割系统。
本发明是这样实现的，提出一种独立成分分析多形状先验水平集方法，实现的技术思路是：首先将形状先验用独立成分分析映射到低维空间，然后在低维空间用Parzen窗估计形状先验的分布，进而设计形状约束项，构造能量函数进行优化求解。所述独立成分分析多形状先验水平集方法包括：输入待分割图像和形状先验；曲线初始化；形状先验对齐；将对齐的形状先验用水平集函数编码；构成形状先验矩阵；用独立成分分析降维；将当前水平集函数投影到低维空间；估计形状先验的概率分布，构造形状驱动能量项，和数据驱动能量项结合，构成能量函数；最小化能量函数，驱动曲线演化，得到分割结果。
进一步，所述独立成分分析多形状先验水平集方法包括以下步骤：
(1)输入大小为M×N的待分割图像I和n个形状先验；
(2)人工或者自动地初始化演化曲线，并用水平集函数编码为φ0；
(3)用基于矩的方法将n个形状先验对齐；
(4)将对齐的形状先验用水平集函数编码，得到
(5)将形状先验的水平集函数按列展开成列向量构成形状先验矩阵
(6)用独立成分分析对形状先验矩阵D进行降维：对形状先验矩阵预处理；对预处理后的形状先验矩阵，用FastICA算法进行ICA降维；
(7)将当前水平集函数φt按列展开成向量，用(6)中得到的投影矩阵将其映射到低维空间；
(8)在低维空间估计形状先验的概率分布，构造形状驱动能量项，并和基于区域的数据驱动能量项结合，构成总的能量函数；
(9)最小化能量函数，驱动曲线演化，得到分割结果。
进一步，所述对形状先验矩阵预处理的方法包括：
第一步，数据中心化：形状先验矩阵减去均值，得到中心化后的形状先验矩阵
第二步，对中心化的形状先验矩阵D1进行主成分分析：首先求D1的协方差矩阵，C＝D1D1T，其中D1T为矩阵D1的转置，对协方差矩阵C进行特征值分解，C＝UΛUT，其中为C的特征值构成的对角矩阵，U＝[u1,...,uM×N]为C的特征向量；
第三步，对中心化的形状先验矩阵D1进行白化处理：计算白化矩阵W＝inv(sqrt(Λ))UT，得到白化后的矩阵满足每个特征维度上的方差为1，并且不同特征维度之间不相关。
进一步，所述对预处理后的形状先验矩阵，用FastICA算法进行ICA处理的步骤如下：
第一步，设定要提取的独立成分的数量，降维之后的特征数k；
第二步，假设映射矩阵为其中，列向量fi∈RM×N,i＝1,...,k，则降维后的形状先验特征向量为构成矩阵Ψ＝[ψ1,...,ψn]；
第三步，用峰度度量降维后形状先验特征向量ψi,i＝1,...,n的非高斯性；
第四步，用固定点迭代算法求解使得峰度最大的方向{f1T,...,fnT}；
第五步，得到映射矩阵F和降维之后的形状先验矩阵
进一步，所述将向量化的当前水平集函数φt降维的步骤包括：
第一步，数据中心化(去均值)，
第二步，将投影到低维空间，
进一步，所述构造能量函数的步骤包括：
第一步，在低维空间用Parzen窗估计形状先验的概率密度函数其中Kσ是方差为σ2高斯核函数，H为Heaviside函数；
第二步，用估计的概率密度的函数构造形状驱动能量项：
第三步，将形状驱动项和基于区域的数据项结合到一起构成总的能量函数：
Etotal＝Edata+λEshape；
其中，λ是相对于数据项，形状约束项的权值。
进一步，所述最小化能量函数，驱动曲线演化的步骤包括：
第一步，形状驱动能量项对应的演化方程为：
第二步，结合形状驱动项和数据驱动项，水平集函数更新方程为：
第三步，迭代更新水平集函数φ，直至收敛得到最后的分割结果φ*。
本发明首先用独立成分分析将形状先验映射到低维空间，再在低维空间中用Parzen窗估计形状先验的统计分布，然后设计形状驱动项，进而构造能量函数进行优化求解。由于形状先验的高维特征存在冗余，并且满足非高斯性，而独立成分分析通过分析高阶统计量可以消除观察信号中高阶统计关联，找出相互独立的隐含信息成分，所以使得本发明对形状先验特征降维的结果更加全面揭示形状先验数据间的本质结构。另外针对未知分布，Parzen窗是一种有效的估计概率密度函数的非参数方法，本发明在用独立成分分析降维后的低维空间用Parzen窗估计形状先验的概率分布，使得本发明估计的形状先验的分布更为准确，从而能够更有效地约束水平集函数的演化，使分割结果更准确。
附图说明
图1是本发明实施例提供的独立成分分析多形状先验水平集方法流程图。
图2是本发明实施例提供的独立成分分析多形状先验水平集方法实现流程图。
图3是本发明实施例提供的在“海星”自然图像中，本发明的分割结果和Cremers等人的多先验水平集分割方法分割结果的对比示意图。
图4是本发明实施例提供的图3对应的形状先验示意图。
图5是本发明实施例提供的在左心室核磁共振医学图像中，本发明的分割结果和Cremers等人的多先验水平集分割方法分割结果的对比示意图。
图6是本发明实施例提供的图5对应的形状先验示意图。
图7是本发明实施例提供的在膝盖骨核磁共振医学图像中，本发明的分割结果和Cremers等人的多先验水平集分割方法分割结果的对比示意图。
图8是本发明实施例提供的图7对应的形状先验示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明能够消除形状先验高维冗余的特征，进而能够更为准确地统计形状先验的分布，形成更有效的形状约束，最后可以得到比较准确的分割结果。
下面结合附图对本发明的应用原理作详细的描述。
如图1所示，本发明实施例提供的独立成分分析多形状先验水平集方法包括以下步骤：
S101：输入待分割图像和形状先验；
S102：曲线初始化；
S103：形状先验对齐；
S104：将对齐的形状先验用水平集函数编码；
S105：构成形状先验矩阵；
S106：用独立成分分析降维；
S107：将当前水平集函数投影到低维空间；
S108：估计形状先验的概率分布，构造形状驱动能量项，和数据驱动能量项结合，构成能量函数；
S109：最小化能量函数，驱动曲线演化，得到分割结果。
下面结合附图对本发明的应用原理作进一步的描述。
如图2所示，本发明实施例提供的独立成分分析多形状先验水平集方法包括以下步骤：
步骤一，输入图像。大小为M×N的待分割图像I和n个形状先验；
步骤二，人工或者自动地初始化演化曲线，并用水平集函数编码为φ0。
步骤三，用基于矩的方法将n个形状先验对齐；
步骤四，将对齐的形状先验用水平集函数编码，得到
步骤五，将形状先验的水平集函数按列展开成列向量可构成矩阵
步骤六，用独立成分分析(ICA)对形状先验矩阵D进行降维；
对形状先验矩阵预处理：
第一步，数据中心化(去均值)：将形状先验矩阵减去均值，即
第二步，对中心化的形状先验矩阵D1进行主成分分析(PCA)：首先求D1的协方差矩阵，C＝D1D1T，其中D1T为矩阵D1的转置，然后对协方差矩阵C进行特征值分解，C＝UΛUT，其中为C的特征值构成的对角矩阵，U＝[u1,...,uM×N]为C的特征向量；
第三步，对中心化的形状先验矩阵D1进行白化处理：计算白化矩阵W＝inv(sqrt(Λ))UT，得到白化后的矩阵满足每个特征维度上的方差为1，并且不同特征维度之间不相关。
对预处理后的形状先验矩阵，用FastICA算法进行ICA处理：
第一步，设定要提取的独立成分的数量，即降维之后的特征数k；
第二步，假设映射矩阵为其中，列向量fi∈RM×N,i＝1,...,k，则降维后的形状先验特征向量为构成矩阵Ψ＝[ψ1,...,ψn]；
第三步，用峰度度量降维后形状先验特征向量ψi,i＝1,...,n的非高斯性；
第四步，用固定点迭代算法求解使得峰度最大的方向{f1T,...,fnT}；
第五步，得到映射矩阵F和降维之后的形状先验矩阵
步骤七，将当前水平集函数φi按列展开成向量并降维；
第一步，数据中心化(去均值)，
第二步，将投影到低维空间，
步骤八，在低维空间估计形状先验的概率分布，构造形状驱动能量项，和基于区域的数据驱动能量项结合，构成能量函数；
第一步，在低维空间用Parzen窗估计形状先验的概率密度函数其中Kσ是方差为σ2高斯核函数，H为Heaviside函数；
第二步，用估计的概率密度的函数构造形状驱动能量项：
第三步，将形状驱动项和基于区域的数据项结合到一起构成总的能量函数：
Etotal＝Edata+λEshape；
其中，λ是相对于数据项，形状约束项的权值。
步骤九，最小化能量函数，驱动曲线演化，得到分割结果。
第一步，形状驱动能量项对应的演化方程为：
第二步，结合形状驱动项和数据驱动项，水平集函数更新方程为：
第三步，迭代更新水平集函数φ，直至收敛得到最后的分割结果φ*。
下面结合仿真实验对本发明的应用效果作详细的描述。
1.仿真条件
本发明是在中央处理器为Intel(R)Corei7-47903.60GHZ、内存16G、WINDOWS7操作系统上，运用MATLAB软件进行的仿真。
2.仿真内容
本发明在自然图像和医学图像上进行了图像分割仿真实验，选取了1幅自然图像和2幅医学图像，并选取Cremers等人的多先验水平集分割方法进行对比，来说明本发明的分割效果。图3是“海星”图像的分割，这幅图像的部分背景和待分割的目标颜色相似，图4是其对应的形状先验。图5是“膝盖骨”核磁共振(MRI)图像的分割，这幅图像待分割的目标和周围组织颜色分布一致，图6是其对应的形状先验。图7是“左心室”核磁共振(MRI)图像的分割，这幅图像待分割的目标和其他器官颜色相似，图8是其对应的形状先验。
3.仿真效果分析：
本发明可以适用于自然图像和医学图像的目标分割。在图像的背景较复杂，和待分割的目标颜色分布相似，用数据项无法正确分割的情况下，本发明可以取得较好的分割结果。说明独立成分分析降维得到的结果可以揭示数据间的本质结构，在低维空降获得的形状先验的分布更为准确，从而能形成有效的形状先验约束项，最后得到较好的分割结果。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。