CN104835168A - 基于全局凸优化变分模型的快速多相图像分割方法 - Google Patents

Abstract

一种基于全局凸优化变分模型的快速多相图像分割方法，主要包括多标记特征函数的定义、非凸能量泛函的构、建能量泛函的转换、全局凸优化的能量泛函构建和基于能量最小化的对偶方法求解过程。在基于区域竞争模型分割模型和多标记特征函数定义的基础上，通过构造非凸的能量泛函和能量泛函的凸化表示，并用对偶方法计算能量泛函的最小化解。本发明所提出的分割方法既能解决非凸目标泛函的局部极优问题，使得分割结果与初始条件无关，又能极大地提高算法的计算效率。

Description

基于全局凸优化变分模型的快速多相图像分割方法

技术领域

本发明涉及的是一种图像处理技术领域图像分割的方法，具体是一种基于全局凸优化变分模型(Global Convex Variational，GCV)的快速多相图像分割方法。

技术背景

图像分割是图像特征提取和分类的重要环节，图像分割的目的就是将图像中的灰度同质区域分离开，并通过各个同质区域的边界来表达。基于变分偏微分方程的图像分割方法，由于具有将模型的初始估计和图像数据先验知识统一于特征提取过程中，并利用分割过程中获得的先验知识来指导分割过程等优势，因此成为近年来图像分割的研究热点。传统的两相Chan-Vese分割模型无法分割多区域同质目标，其非凸能量泛函使得目标函数陷入局部最优，导致分割结果依赖于演化曲线的初始化，而且还使得一些快速数值计算方法无法应用。基于变分水平集的多相图像分割方法，采用差分格式求解复杂的目标函数导致计算效率非常低，极大地限制了该算法的实际应用价值。因此，如何建立有效分割模型的能量泛函，以消除非凸能量泛函局部极优问题使得分割结果依赖于初始条件的影响，是多相图像分割领域的研究热点之一。

经对现有技术文献的检索发现，基于变分法理论的多相水平集方法最初通过集成多种信息来构建分割模型(Zhao H.K.“A variational level set approach to multiphase motion”(多相运动的变分水平集方法)，Journal of Computational Physics(1996)7:79-195)。多相分段恒定的变分水平集分割模型(Vese L.,Chan T.“A multiphase level set framework for imagesegmentation using the Mumford and Shah model”(基于Mumford-Shah模型的多相水平集图像分割框架)，International Journal of Computer Vision)(2002)50(3):271-293)，进一步简化了Mumford-Shah模型，用N个水平集函数表达2^N个相，避免多个水平集函数的重叠。方江雄等提出的多分辨率多区域变分水平集图像分割方法(专利号：CN102044077B)，通过多分辨率技术执行分割图像的演化曲线来解决初始化水平集能量函数陷入局部能量最小值，降低了噪声的干扰、减小了搜索的空间。但是，这些多相图像分割方法所建立的能量泛函均是非凸函数，导致局部极优问题的产生，而且还使得现有的快速算法无法直接应用。

发明内容

本发明的目的是，通过提出一种基于全局凸优化变分模型的快速多相图像分割方法，来解决非凸目标泛函的局部极优问题，使得分割结果与初始条件无关，并能用快速的数值计算方法求解。

本发明的技术方案：在基于区域竞争模型分割模型和多标记特征函数定义的基础上，通过构造非凸的能量泛函和能量泛函的凸化表示，并用对偶方法计算能量泛函的最小化问题，从而避免能量泛函的局部极优问题。具体步骤如下：

步骤1：多标记特征函数的定义。在不相交的图像子域Ω_i中(整个图像域为Ω)的点x，多标记特征函数u_i(x)∈[0,1]，其中i＝1,……,N，其表达式如下：

$u_i\left(x\right)=1_{{\mathrm{\Ω}}_i}=\left\{\begin{array}{cc}1& x\∈{\mathrm{\Ω}}_i\\ 0& x\∉{\mathrm{\Ω}}_i\end{array}\right.,i=1,......,N$

多标记约束项满足

步骤2：非凸能量泛函的构建。本发明采用了Fang Jiangxiong(方江雄)等在OpticalEngineering上Statistical approaches to automatic level set image segmentation with multipleregions中定义的能量泛函。假设给定的图像子域Ω中点坐标为x，用N-1个水平集函数φ_i表达N个区域(i＝1,…,N-1)，其能量泛函表达式如下：

$\underset{(r,\mathrm{\φ})}{\min }E(r,\mathrm{\φ})=\underset{(c,\mathrm{\φ})}{\min }\left\{\begin{array}{c}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}|\▿H\left({\mathrm{\φ}}_i\left(x\right)\right)|\mathrm{dx}+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_i)}^{2}H\left({\mathrm{\φ}}_i\left(x\right)\right)\underset{k=0}{\overset{k=i-1}{\mathrm{\Π}}}(1-H\left({\mathrm{\φ}}_i\left(x\right)\right))\mathrm{dx}\\ +{\mathrm{\λ}}_N{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_N)}^2\underset{k=1}{\overset{N-1}{\mathrm{\Π}}}(1-H\left({\mathrm{\φ}}_i\left(x\right)\right))\mathrm{dx}\end{array}\right\}$

其中r＝(r₁,…,r_N),I为图像的灰度均值，表示图像梯度，φ＝(φ₁,…,φ_N-1)，H(x)是Heaviside函数，H(φ₀)≡0，δ(x)为平滑函数H(x)的导数，λ₁＞0,i＝1,…,N。根据欧拉拉格朗日等式最小化能量泛函，其演化方程如下：

其中t时间变量，假设H_i＝H(φ_i(x))，变量Φ_i(φ_i)和像素均值c_i表达式如下：

$\left\{\begin{array}{c}{\mathrm{\Φ}}_i\left({\mathrm{\φ}}_i\right)={\mathrm{\λ}}_{i+1}{(I-c_{i+1})}^2{H}_{i+1}+{\mathrm{\λ}}_{i+2}{(I-c_{i+2})}^2{\hat{H}}_{i+1}{H}_{i+2}+...\\ +{\mathrm{\λ}}_{N-1}{(I-c_{N-1})}^2{\hat{H}}_{i+1}...{\hat{H}}_{N-2}{H}_{N-1}+{\mathrm{\λ}}_N{(I-c_N)}^2{\hat{H}}_{i+1}...{\hat{H}}_{N-2}{\hat{H}}_{N-1}\\ c_i=\frac{{\∫}_{\mathrm{\Ω}}I\·H_i\underset{k=0}{\overset{k=i-1}{\mathrm{\Π}}}{\hat{H}}_k\mathrm{dx}}{{\∫}_{\mathrm{\Ω}}{H}_i\underset{k=0}{\overset{k=i-1}{\mathrm{\Π}}}{\hat{H}}_k\mathrm{dx}}\end{array}\right.$

为了找到全局最小化分割模型，通过去掉函数δ(x)，其分割模型的能量泛函如下：

$\underset{(c,\mathrm{\φ})}{\min }E(c,\mathrm{\φ})=\underset{(c,\mathrm{\φ})}{\min }\left\{\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}\right|\▿{\mathrm{\φ}}_i\left(x\right)|\mathrm{dx}+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x){\mathrm{\φ}}_i\left(x\right)\mathrm{dx}\}$

其中c＝(c₁,…,c_N)。

步骤3：能量泛函的转换。为求解步骤2中的非凸能量泛函，多标记特征函数重定义如下：

其中λ₀＞0。通过定义变量u_i,i＝1,…,N-1，本发明采用凸松弛方法来解决非凸问题，假设其凸集□定义如下：

用图像变量来替代标记函数，其能量泛函表达式可重写为：

$\underset{U=(u_1,...,u_{N-1})\∈{\left[\mathrm{0,1}\right]}^{N-1}}{\min }(c,U)=\underset{(c,U\∈{\left[\mathrm{0,1}\right]}^{N-1})}{\min }\left\{\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}\right|\▿u_i\left(x\right)|\mathrm{dx}+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}(c_i,{\mathrm{\Φ}}_i,x)u_i\left(x\right)\mathrm{dx}\}$

其中U＝{u₁,…,u_N-1}。该能量泛函由数据项和规则项构成，规则项中用全变分范数表示，函数g(x)＝1/(1+x²)，能量泛函表示可改写为：

$\underset{U\∈\left[\mathrm{0,1}\right]}{\min }\{\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{TV}}_g\left(u_i\left(x\right)\right)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)u_i\left(x\right)\mathrm{dx}\}$

步骤4：全局凸优化的能量泛函构建。针对任意c_i∈R⁺,i＝1,…,N-1，函数g(x)∈[0,1，]步骤3中非凸约束能量泛函转化为凸优化非约束最小化能量泛函如下：

$\underset{U\∈{\left[\mathrm{0,1}\right]}^{N-1}}{\min }\left\{\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}\right|\▿u_i\left(x\right)|\mathrm{dx}+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_i)}^2{u}_i\underset{k=0}{\overset{k=i-1}{\mathrm{\Π}}}{\hat{u}}_k+{\mathrm{\λ}}_N{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_N)}^2\underset{k=1}{\overset{k=N-1}{\mathrm{\Π}}}{\hat{u}}_k\mathrm{dx}\}$

其中常量α＞0远大于λ₂，κ(u_i)＝max{min{2|u_i-1|},1}，惩罚因子α＞1/2||R(x)||_L∞。

步骤5：基于能量最小化的对偶方法求解过程。为最小化能量泛函E(U,c)，本发明用对偶方法增加辅助变量来解决最小化问题。针对每相凸问题变分问题，通过用图像变量进行参数交替，其表达式如下：

$\underset{u_i\∈\left[\mathrm{0,1}\right]}{\min }{E}^{\mathrm{Cov}}(u_i,c_i)=\underset{u_i\∈\left[\mathrm{0,1}\right]}{\min }\{{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)u_i\left(x\right)\mathrm{dx}+\mathrm{\α\κ}\left(v_i\right)+{\mathrm{TV}}_g\left(u_i\right)\}$

各子区域Ω_i像素均值在迭代过程中更新方程式如下：

其中 ${\hat{u}}_i=1-u_i,i=1,...,N-1.$

依据Chan等在SIAM Journal on Applied Mathematics上提出的Algorithms for findingglobal minimizers of image segmentation and denoising models中方法，通过在能量泛函中增加对偶变量(u_i,v_i)，其能量泛函可改写如下：

$\underset{u_i\∈\left[\mathrm{0,1}\right]}{\min }{E}^{\mathrm{Cov}}(u_i,c_i)=\underset{u_i\∈\left[\mathrm{0,1}\right]}{\min }\{{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)v_i\left(x\right)\mathrm{dx}+\mathrm{\α\κ}\left(v_i\right)+\frac{1}{2\mathrm{\θ}}{\left|\right|u_i-v_i\left|\right|}^2+{\mathrm{TV}}_g\left(u_i\right)\}$

其中θ是非常小值。因为能量泛函是凸函数，最小化能量泛函就能得到全局最小值。为解决u_i和v_i的凸优化问题，采用两步迭代算法计算最小化问题。

第一步：当v_i固定，求解v_i，其表达可写成：

$\underset{u_i\left(x\right),v_i\left(x\right)}{\min }\{{\mathrm{TV}}_g\left(u_i\right)+\frac{1}{2\mathrm{\θ}}{\left|\right|u_i-v_i\left|\right|}^2\}$

假定u_i(x)＝v_i(x)-θ_idivp_i，向量p_i满足等式通过固定点方法

$p_i^{n+1}=\frac{p_i^n+\mathrm{\δ}t\▿(\mathrm{div}\left(p_i^n\right)-v_i/\mathrm{\θ})}{1+\frac{\mathrm{\δt}}{g\left(x\right)}|\▿(\mathrm{div}\left(p_i^n\right)-v_i/\mathrm{\θ})|}$

我们通过条件来终止固定点迭代问题。

第二步：求解v_i(x)，与v_i(x)能量泛函表达式如下：

$\underset{v_i\left(x\right)}{\min }\{{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)v_i\left(x\right)\mathrm{dx}+\mathrm{\α\κ}\left(v_i\right)+\frac{{\mathrm{\λ}}_1}{2\mathrm{\θ}}{\left|\right|u_i-v_i\left|\right|}^2\}$

其中v_i(x)＝min{max{u_i(x)-θp_i(x),0},1}。

本发明通过在基于区域竞争模型分割模型基础上，构造非凸的能量泛函和能量泛函的凸化表示，并用对偶方法计算能量泛函的最小化解。本发明所提出的方法既能解决非凸目标泛函的局部极优问题，使得分割结果与初始条件无关，又能极大地提高算法的计算效率。

附图说明

图1表示本发明实施例中基于全局凸优化变分模型的快速多相图像分割方法流程图；

图2为在不同的初始位置条件下，基于全局凸优化变分模型的快速方法分割医学图像结果。其中：第一行表示分割图像的初始轮廓图；第二行和第三行分别表示在不同初始条件下演化曲线的最终停止位置和分割结果；

图3比较了本发明所提出的GCV分割模型和多区域竞争分割(Multi-region CompetitionSegmentation，MCS)模型分割遥感图像结果；

其中：图3(a)和图3(e)分别显示了MCS模型和GCV模型分割遥感图像的初始轮廓；图3(b)和图3(c)分别显示了MCS模型分割结果的两目标区域；图3(f)和图3(g)分别显示了GCV模型分割结果的两目标区域；图3(d)和图3(h)分别显示了MCS模型和GCV模型分割遥感图像结果。

具体实施方式

在基于区域竞争模型分割模型和多标记特征函数定义的基础上，通过构造非凸的能量泛函和能量泛函的凸化表示，并用对偶方法求解能量泛函的最小解。本发明所提出的分割方法既能解决非凸目标泛函的局部极优问题，使得分割结果与初始条件无关，又能极大地提高分割算法的计算效率。本发明具体实施步骤包括如下：

(1)输入分割图像，设置初始化参数：分割区域数目N，两个权重系数附加θ，最大迭代次数；

(2)各区域参量初始化：计算各区域的像素均值c_i，特征函数u_i和v(u_i)，初始化其中i＝1,…N-1；

(3)运行如下的递归过程当满足条件或循环次数不大于最大迭代次数，重复如下操作：

①通过固定点方法所求解的公式，计算各向量p_i,i＝1,…N-1；

②计算各区域像素均值c_i,i＝1,…N；

③根据公式v(u_i)＝min{max{u_i(x)-θ_ip_i(x),0},1}，更新参数v(u_i),i＝1,…N-1。

(4)输出各分割区域的图像和分割结果。

图2显示了在不同的初始位置条件下，基于全局凸优化变分模型的快速算法分割医学图像结果。实验中，分割区域数为3，试验中区域标记参数λ_i＝1，i＝1,…,N。图中第一行分别显示了初始化轮廓形状分别为两个圆形、矩形、直线和两个三角形。第二行和第三行分别显示了基于GCV模型的多相分割方法分割同一幅医学图像后，演化曲线最终停止位置和分割结果。

图3比较了本发明所提出的GCV分割模型和多区域竞争分割(MCS)模型分割遥感图像结果。实验中，分割区域数为3，区域标记参数λ₀＝0.5。在两组图像试验中，图像初始化轮廓位置相同，图(a)和图(e)显示了演化曲线最终停靠位置轮廓。图(b)和图(c)显示MCS模型对应的两目标区域。图(f)和图(g)显示GCV模型对应的两目标区域。图(d)和图(h)显示了MCS模型和GCV模型分割的最终结果。从分割的结果看，基于GCV模型的分割方法比基于MCS模型的分割方法有更好的效果。

Claims (4)

1.一种基于全局凸优化变分模型的快速多相图像分割方法，其能量泛函构建过程如下：

假设给定的图像子域Ω中点坐标为x，用N-1个水平集函数φ_i表达N个区域(i＝1,…,N-1)，其能量泛函表达式如下：

$\underset{(r,\mathrm{\φ})}{\min }E(r,\mathrm{\φ})=\underset{(c,\mathrm{\φ})}{\min }\left\{\begin{array}{c}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}|\▿H\left({\mathrm{\φ}}_i\left(x\right)\right)|\mathrm{dx}+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_i)}^{2}H\left({\mathrm{\φ}}_i\left(x\right)\right)\underset{k=0}{\overset{k=i-1}{\mathrm{\Π}}}(1-H\left({\mathrm{\φ}}_i\left(x\right)\right))\mathrm{dx}\\ +{\mathrm{\λ}}_N{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_N)}^2\underset{k=1}{\overset{N-1}{\mathrm{\Π}}}(1-H\left({\mathrm{\φ}}_i\left(x\right)\right))\mathrm{dx}\end{array}\right\}$

其中r＝(r₁,…,r_N),I为图像的灰度均值，表示图像梯度，φ＝(φ₁,…,φ_N-1)，H(x)是Heaviside函数，H(φ₀)≡0，δ(x)为平滑函数H(x)的导数，λ₁＞0,i＝1,…,N，其分割方法的具体步骤如下：

步骤1：输入分割图像，定义多标记特征函数；在不相交的图像子域Ω_i中，给每个区域定义一个多标记特征函数u_i(x)∈[0,1](i＝1,……,N-1)；

步骤2：分割模型能量泛函的构建，根据图像概率密度函数分布规律，用N-1个水平集函数φ_i表达N个区域(i＝1,…,N-1)；

步骤3：能量泛函的转换，重新定义多标记特征函数，将规则项中用全变分范数 ${\mathrm{TV}}_g\left(u_i\right)={\∫}_{\mathrm{\Ω}}g\left(\right|\▿I\left(x\right)\left|\right)\left|{\▿u}_i\right|\mathrm{dx}$ 表示，来简化能量泛函形式；

步骤4：全局凸优化能量泛函的构建，将非凸约束能量泛函转化为凸优化非约束最小化能量泛函；

步骤5：基于能量最小化的对偶方法求解过程，通过增加辅助变量v_i(x)，用对偶方法来解决能量泛函的最小化问题。

2.根据权利要求1所述的基于全局凸优化变分模型的快速多相图像分割方法，其特征在于：能量泛函的转换，多标记特征函数重定义如下：

其中λ₀＞0，通过已定义的变量u_i,i＝1,…,N-1，其中采用凸松弛方法来解决非凸问题，其凸集□定义如下：

用图像变量来替代标记函数，其能量泛函表达式可重写为：

$\underset{U=(u_1,...,u_{N-1})\∈{[0,i]}^{N-1}}{\min E(c,U)}=\underset{(c,U\∈{\left[\mathrm{0,1}\right]}^{N-1})}{\min }\left\{\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}\right|{\▿u}_i\left(x\right)|\mathrm{dx}+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)u_i\left(x\right)\mathrm{dx}\}$

其中U＝{u₁,…,u_N-1}，该能量泛函由数据项和规则项构成，规则项中用全变分范数表示，函数g(x)＝1/(1+x²)，能量泛函表示可改写为：

$\underset{U\∈\left[\mathrm{0,1}\right]}{\min }\{\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{TV}}_g\left(u_i\left(x\right)\right)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)u_i\left(x\right)\mathrm{dx}\}.$

3.根据权利要求1所述的基于全局凸优化变分模型的快速多相图像分割方法，其特征在于：全局凸优化的能量泛函构建，对任意c_i∈R⁺,i＝1,…,N-1，函数g(x)∈[0,1]，步骤3中非凸约束能量泛函转化为凸优化非约束最小化能量泛函如下：

$\underset{U\∈{\left[\mathrm{0,1}\right]}^{N-1}}{\min }\left\{\begin{array}{c}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\∫}_{\mathrm{\Ω}}\left|{\▿u}_i\left(x\right)\right|\mathrm{dx}+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_i)}^2{u}_i\underset{k=0}{\overset{k=i-1}{\mathrm{\Π}}}{\hat{u}}_k\\ +{\mathrm{\λ}}_N{\∫}_{\mathrm{\Ω}}{(I\left(x\right)-c_N)}^2\underset{k=1}{\overset{k=N-1}{\mathrm{\Π}}}{\hat{u}}_k\mathrm{dx}\end{array}\right\}$

其中常量α＞0远大于λ₂，κ(u_i)＝max{min{2|u_i-1|},1}，惩罚因子

4.根据权利要求1所述的基于全局凸优化变分模型的快速多相图像分割方法，其特征在于：基于能量最小化的对偶方法求解过程，通过用对偶方法增加辅助变量来解决最小化问题，其表达式如下：

$\underset{u_i\∈[0,i]}{\min }{E}^{\mathrm{Cov}}(u_i,c_i)=\underset{u_i\∈[0,i]}{\min }\{{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)u_i\left(x\right)\mathrm{dx}+\mathrm{\α\κ}\left(v_i\right)+{\mathrm{TV}}_g\left(u_i\right)\}$

各子区域Ω_i像素均值在迭代过程中更新方程式如下：

其中 ${\hat{u}}_i=1-u_i,i=1,...,N-1,$

依据Chan等在SIAM Journal on Applied Mathematics上提出的Algorithms for findingglobal minimizers of image segmentation and denoising models中方法，通过在能量泛函中增加对偶变量(u_i,v_i)，其能量泛函可改写如下：

$\underset{u_i\∈\left[\mathrm{0,1}\right]}{\min }{E}^{\mathrm{Cov}}(u_i,c_i)=\underset{u_i\∈\left[\mathrm{0,1}\right]}{\min }\{{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)v_i\left(x\right)\mathrm{dx}+\mathrm{\α\κ}\left(v_i\right)+\frac{1}{2\mathrm{\θ}}{\left|\right|u_i-v_i\left|\right|}^2+{\mathrm{TV}}_g\left(u_i\right)\}$

其中θ是非常小值，因为能量泛函是凸函数，最小化能量泛函就能得到全局最小值，为解决u_i和v_i的凸优化问题，采用两步迭代算法计算最小化问题，

第一步：当v_i固定，求解v_i，其表达可写成：

$\underset{u_i\left(x\right),v_i\left(x\right)}{\min }\{{\mathrm{TV}}_g\left(u_i\right)+\frac{1}{2\mathrm{\θ}}{\left|\right|u_i-v_i\left|\right|}^2\}$

假定u_i(x)＝v_i(x)-θ_idivp_i，向量p_i满足等式 $g\left(x\right)\▿({\mathrm{\θ}}_i{\mathrm{div}p}_i-v_i)-\left|({\mathrm{\θ}\mathrm{div}p}_i-v_i)\right|p_i=0,$ 通过固定点方法

$p_i^{n+1}=\frac{p_i^n+\mathrm{\δt}\▿(\mathrm{div}\left(p_i^n\right)-v_i/\mathrm{\θ})}{1+\frac{\mathrm{\δt}}{g\left(x\right)}|\▿(\mathrm{div}\left(p_i^n\right)-v_i/\mathrm{\θ}|}$

我们通过条件 $|p_i^{n+1}\left(x\right)-p_i^n\left(x\right)|\≤\mathrm{\β},(\mathrm{\β}={10}^{-2})$ 来终止固定点迭代问题，

第二步：求解v_i(x)，与v_i(x)能量泛函表达式如下：

$\underset{v_i\left(x\right)}{\min }\{{\∫}_{\mathrm{\Ω}}r(c_i,{\mathrm{\Φ}}_i,x)v_i\left(x\right)\mathrm{dx}+\mathrm{\α\κ}\left(v_i\right)+\frac{{\mathrm{\λ}}_1}{2\mathrm{\θ}}{\left|\right|u_i-v_i\left|\right|}^2\}$

其中v_i(x)＝min{max{u_i(x)-θp_i(x),0},1}。