融合变分和深度先验的非均匀图像分割算法

Heterogeneous Image Segmentation Algorithm Combining Variational and Depth Prior

摘要: 为了尽可能精确地分割出强度不均匀图像更多的细节部分，提出一种融合深度图像先验的变分图像分割模型，利用交替方向乘子法设计相应的数值求解算法。实验结果表明，提出的模型在去噪正则化(regularization by denoising, RED)框架下融合了TV正则项捕获边缘和卷积神经网络(convolutional neural network, CNN)捕获细节的优势，尤其在处理结构丰富和纹理细致的图像时，可以分割出更多的细节，分割结果更精确。同时，提出的方法对于不同的初始轮廓具有很好的鲁棒性。此外，与对比实验中处理非均匀图像分割的方法相比，该模型算法复杂度低，具有快速高效的优势。

Abstract: To segment more detailed parts of a homogeneous image as accurately as possible, a variational image segmentation model with depth image prior is proposed in this paper. The alternating direction method of multipliers (ADMM) is used to design the corresponding numerical algorithm. The results reveal that the proposed model by incorporating the denoising regularization (RED) framework with TV regularization can retain more details in images and obtain more accurate segmentation accuracy. At the same time, the model is robust to different initial contours. In addition, compared with the image segmentation method in the comparison experiment, the algorithm complexity of the proposed method is lower and has higher operational efficiency.

1. 引言

图像分割作为图像处理和计算机视觉领域中的重要任务，其主要目标是将图像划分成若干具有相似特征的区域。然而，现实世界中的图像由于照明变化和成像设备的影响，常常呈现出强度不均匀性。这对传统基于强度均匀性假设的分割方法提出了挑战，使得对具有不同强度的图像进行准确分割成为一项关键而复杂的任务。

在过去的二十年里，变分方法和偏微分方程(partial differential equation, PDE)已成功应用于图像分割领域。这些方法大致可分为两类：基于边缘的方法[1] [2]和基于区域的方法[3] [4]。基于边缘的方法涵盖了通过演化初始轮廓来识别物体边界的方法，尽管这些方法能够有效地分割图像，但对初始轮廓和目标边界敏感。而基于区域的方法则利用区域信息引导轮廓演化。然而，大多数基于区域的模型都假设强度均匀性，因此，这些模型在处理强度不均匀的图像时表现不佳。目前，已经出现了几种基于Retinex理论的方法解决由强度不均匀性带来的图像分割挑战[5]-[8]。2020年，Jin等人[9]在CVB模型[8]中引入TV正则化器描述反射部分，提出了CVBS模型。CVBS模型提高了图像轮廓的分割精度，可以产生相对精确的分割结果。基于Retinex理论的方法在对强度不均匀的图像分割时，能够较好地保持物体的轮廓，但是在分割一些纹理鲜明的图像时，往往不能较精准地将一些细节部分分割出来，在这些区域会出现明显的欠分割和边界模糊等问题。

近年来，随着深度学习的发展与完善，卷积神经网络(Convolutional Neural Network，CNN)可以从图像中学习到更复杂的特征，从而捕获更精确的细节。考虑到这一因素，越来越多的研究注重于将深度图像先验与传统的图像处理方法融合。目前，即插即用(Plug-and-Play, PnP)框架[10]-[13]和去噪正则化(Regularization by Denoising, RED)框架[14]-[16]是两个最常用的融合框架。根据PnP理论，去噪器可以作为基于模型方法的隐式图像先验，用于解决多数逆问题，如图像恢复、图像分割等。然而，在即插即用框架中，隐式正则项缺乏物理解释，需要进行苛刻的参数调整，导致算法复杂度高，运行效率缓慢。为了解决这些问题，Romano等人[14]提出了一个更强大、更灵活的替代框架——去噪正则化框架。该框架使用去噪引擎定义逆问题的正则化，提出一种显式的图像正则化函数。RED可以完全灵活地选择最小化泛函的迭代优化过程，可以结合图像去噪算法，非常有效地处理一般逆问题，并且参数较少，运行速率较高。

本文在RED框架[17]-[20]下把深度图像先验与CVBS模型结合，提出了一种融合深度图像先验的分割模型。主要工作如下：① 深度神经网络可以学习更复杂的特征表示，更好地捕获图像中的语义信息，而传统TV正则边缘保持较好，本文利用RED框架将二者融合，所提模型可以明显地提高图像分割精度；② 使用RED框架，与一些传统图像分割模型相比，算法复杂度低，大大减少了运行时间，提高了运行效率；③ 优化问题对于单个变量是凸的，因此可以使用交替方向乘子法(alternating direction method of multipliers, ADMM)，这种方法确保所提的方法对于不同的初始轮廓具有很好地鲁棒性。

2. 相关工作

本节主要介绍论文相关的理论知识和模型，包括RED框架、基于Retinex理论的CVBS模型。

2.1. RED框架

RED框架是一种将外部去噪方法融入到重建过程中的图像重建框架[17]。其核心概念是利用现有的去噪算法作为正则化项，以帮助解决成像中的逆问题。该框架使用以下表达式作为正则化项：

$R\left(u\right)=\frac{1}{2}{u}^{\text{T}}\left[u-f\left(u\right)\right],$ (1)

其中$f(\cdot )$ 是选择的去噪器，$u$ 是带噪声的图像，并且

$\nabla R\left(u\right)=u-f\left(u\right).$ (2)

RED框架依赖于一种图像自适应拉普拉斯算子的形式，它建立了一个强大的先验，可用于正则化逆问题。有两个重要的性质：① 先验的梯度由式(2)直接给出，避免了对去噪函数$f\left(u\right)$ 的微分。② $R\left(u\right)$ 是一个凸函数，去噪引擎可以以不同的方式插入到迭代方案中。所有这些方案都涉及与数据一致性项相关的子问题和包含与去噪子问题相关的去噪引擎的正则化项。因此，可以使用许多高效的去噪算法来解决逆问题。

2.2. 基于Retinex理论的CVBS模型

由Land和McCann [21] [22]提出的Retinex理论有助于分割强度不均匀图像。它将真实的色彩与反射分离开来，通过将图像表现为反射(固有的物体颜色)与照明(光的相互作用)两部分来消除光照和阴影。数学上，这种关系可以表示为

$i_0\left(x\right)=s\left(x\right)b\left(x\right),$ (3)

利用对数变换，

$I_0=\log \left(i\right),B=\log \left(b\right),S=\log \left(s\right),$ (4)

可以得到

$I_0=B+S,$ (5)

为了在强度不均匀性的情况下提高图像分割的精度，Jin [11]提出了如下CVBS模型：

$\underset{{}_{\begin{array}{c}0\le u\le 1\\ B,\overline{S},c_1,c_2\end{array}}}{\min }E\left(u,B,S,c_1,c_2\right),\text{s}\text{.t}{I}_0=B+S,$ (6)

能量泛函为

$\begin{array}{c}E\left(u,B,S,c_1,c_2\right)={\lambda }_2{\displaystyle {\int }_{\Omega }{\left(c_2-S\right)}^2}\left(1-u\right)\text{d}x+{\lambda }_1{\displaystyle {\int }_{\Omega }{\left(c_1-S\right)}^{2}u\text{d}x}+{\displaystyle {\int }_{\Omega }\left|Du\right|}\\ +\alpha {\displaystyle {\int }_{\Omega }{\left|\nabla B\right|}^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }\left|DS\right|},\end{array}$ (7)

其中${\lambda }_1,{\lambda }_2,\alpha ,\beta$ 为正参数，$I_0$ 是观测到的图像，由反射率部分$S$ 和偏置部分$B$ 组成，$u$ 为松弛的Heaviside函数，满足$u\in \left[0,1\right]$ ，$c_1,c_2$ 为两个常数，分别用于近似分割曲线内外的图像强度。前三项是作用于$S$ 的松弛CV模型[3]，第四项保证了整个域上偏置场的平滑性，最后一项加入了TV正则化器表示反射分量s。

CVBS模型可以很好地对简单的非均匀图像进行分割，但在处理一些复杂的图像时，往往不能较精确地分割出物体的某些特征。例如，在图1中，图1(a)和图1(f)分别为原始图像和标准分割结果；图1(b)为CVBS模型的分割结果，从图1(b)可以看出CVBS模型可以大致分割出鹰的轮廓，但是在鹰右翼处的内部(绿色方框区域)，其分割结果存在着欠分割的现象。另外，从图1(c)可以看到，CVBS模型并不能较为准确地分割出鹰每个翅尖结构。

(a)

(b) (c)

(DSC = 0.9584, HD =1.2511)

(d) (e)

(DSC = 0.9944, HD = 0.4045)

(f) (g)

Figure 1. Segmentation results of the two models. (a) Original image; (b) CVBS model; (c) (b)’s enlarged picture; (d) The proposed model; (e) (d)’s enlarged picture; (f) Ground truth; (g) (f)’s enlarged picture

图1. 两种模型的分割结果。(a) 原始图像；(b) CVBS模型；(c) (b)放大的图；(d) 本文模型；(e) (d)放大的图；(f) 标准分割结果；(g) (f)放大的图

3. 本文模型及其快速算法

本节主要介绍本文所提出的模型，并给出模型的算法。

3.1. 本文模型

为了解决CVBS模型不能精确地分割出非均匀图像的某些细节部分出现的欠分割的情况，本文在RED框架下将深度图像先验融合到CVBS模型中，提出了一种融合深度图像先验的分割模型。更准确地说，考虑以下变分问题，

$\underset{{}_{\begin{array}{c}0\le u\le 1\\ B,\overline{S},c_1,c_2\end{array}}}{\min }E\left(u,B,S,c_1,c_2\right),\text{s}\text{.t}{I}_0=B+S,$ (8)

能量泛函为

$\begin{array}{c}E\left(u,B,S,c_1,c_2\right)={\displaystyle {\int }_{\Omega }\left|Du\right|}+\alpha {\displaystyle {\int }_{\Omega }{\left|\nabla B\right|}^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }\left|DS\right|}+\frac{{\lambda }_3}{2}u\left(u-f\left(u\right)\right)\\ +{\lambda }_1{\displaystyle {\int }_{\Omega }{\left(c_1-S\right)}^{2}u\text{d}x}+{\lambda }_2{\displaystyle {\int }_{\Omega }{\left(c_2-S\right)}^2}\left(1-u\right)\text{d}x,\end{array}$ (9)

其中，${\lambda }_1,{\lambda }_2,{\lambda }_3,\alpha ,\beta$ 是正参数，$f(\cdot )$ 是去噪器，本文所使用的是去噪卷积神经网络FFDNet [23]。

相较于CVBS模型，本文方法有效地整合了TV正则化器保持边缘的优势以及CNN捕捉细节的优势，尤其在处理结构丰富和纹理细致的图像时，分割效果较好。如图1(d)所示，可以看出，本文方法将鹰翅膀内部完整地分割出来，并且很好地保持了每一个翅尖原有的结构。同时，定量指标(DSC值和HD值)也说明了本文方法对非均匀图像的分割效果较好。

3.2. 本文算法

本节利用RED框架快速灵活的特点，结合ADMM算法[24] [25]和变量分裂技术，提出了模型(9)的数值求解算法。通过引入辅助变量$d_1,d_2,W$ ，模型(9)等价于以下约束优化问题，

$\begin{array}{l}\underset{{}_{u,B,S}}{\min }\left\{{\displaystyle {\int }_{\Omega }\left|\left|d_1\right|\right|\text{d}x}+\alpha {\displaystyle {\int }_{\Omega }{\left|\nabla B\right|}^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }\left|d_2\right|\text{d}x}+{\lambda }_1{\displaystyle {\int }_{\Omega }{\left(c_1-S\right)}^{2}W\text{d}x}+\frac{{\lambda }_3}{2}u\left(u-f\left(u\right)\right)\right\}\\ +{\lambda }_2{\displaystyle {\int }_{\Omega }{\left(c_2-S\right)}^2}\left(1-W\right)\text{d}x,\\ \text{s}\text{.t}.I_0=B+S,\nabla u=d_1,\nabla S=d_2,u=W.\end{array}$ (10)

利用拉格朗日乘子${\left[\theta ,{\theta }_1,{\theta }_2,{\theta }_3\right]}^{\text{T}}$ ，则式(10)的增广拉格朗日函数为

$\begin{array}{l}ℒ\left(u,S,B,d_1,d_2,\theta ,{\theta }_1,{\theta }_2,{\theta }_3\right)\\ ={\displaystyle {\int }_{\Omega }\left|d_1\right|\text{d}x}+\beta {\displaystyle {\int }_{\Omega }\left|d_{\text{2}}\right|\text{d}x}+\alpha {\displaystyle {\int }_{\Omega }{\left|\nabla B\right|}^2\text{d}x}+{\lambda }_1{\displaystyle {\int }_{\Omega }{\left(c_1-S\right)}^{2}W\text{d}x}+\frac{{\lambda }_3}{2}u\left(u-f\left(u\right)\right)\\ +{\lambda }_2{\displaystyle {\int }_{\Omega }{\left(c_2-S\right)}^2}\left(1-W\right)\text{d}x+\frac{\rho }{2}{\displaystyle {\int }_{\Omega }{\left|I_0-B-S\right|}^2\text{d}x}+\langle \theta ,I_0-B-S\rangle \\ +\frac{{\rho }_1}{2}{\displaystyle {\int }_{\Omega }{\left|d_1-\nabla u\right|}^2\text{d}x}+\langle {\theta }_1,\nabla u-d_1\rangle +\frac{{\rho }_2}{2}{\displaystyle {\int }_{\Omega }{\left|d_2-\nabla S\right|}^2\text{d}x}+\langle {\theta }_2,\nabla S-d_2\rangle \\ +\frac{{\rho }_3}{2}{\displaystyle {\int }_{\Omega }{\left|W-u\right|}^2\text{d}x}+\langle {\theta }_3,W-u\rangle ,\end{array}$ (11)

式中$\rho ,{\rho }_1,{\rho }_2,{\rho }_3$ 为惩罚参数。

给定第k步的迭代，任意中间解$u^k,W^k,S^k,B^k,d_1^k,d_2^k,{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k$ ，用交替极小化法更新问题(11)解的过程如下：

① 固定$u^k,S^k,{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k$ ，更新$d_1^{k+1},d_2^{k+1},B^{k+1}$ 。

$\left(d_1^{k+1},d_2^{k+1},B^{k+1}\right)=\underset{d_1,d_2,B}{\arg \min }L\left(u^k,W^k,S^k,{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k\right)$ (12)

可得

$d_1^{k+1}=\mathrm{shrink}\left(\nabla u^k+\frac{{\theta }_1^k}{{\rho }_1},{\rho }_1\right)=\frac{\nabla u^k+\frac{{\theta }_1^k}{{\rho }_1}}{\left|\nabla u^k+\frac{{\theta }_1^k}{{\rho }_1}\right|}\max \left\{\left|\nabla u^k+\frac{{\theta }_1^k}{{\rho }_1}\right|-\frac{1}{{\rho }_1},0\right\},$ (13)

$d_2^{k+1}=\mathrm{shrink}\left(\nabla S^k+\frac{{\theta }_2^k}{{\rho }_2},{\rho }_2\right)=\frac{\nabla S^k+\frac{{\theta }_2^k}{{\rho }_2}}{\left|\nabla S^k+\frac{{\theta }_2^k}{{\rho }_2}\right|}\max \left\{\left|\nabla S^k+\frac{{\theta }_2^k}{{\rho }_2}\right|-\frac{\beta }{{\rho }_2},0\right\}.$ (14)

接着通过式(15)求解$B^{k+1}$ ，

$\begin{array}{l}{B}^{k+1}=\underset{B}{\arg \min }L\left(u^k,W^k,S^k,d_1^{k+1},d_2^{k+1},{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k\right)\\ =\underset{B}{\arg \min }\alpha {\displaystyle {\int }_{\Omega }{\left|\nabla B\right|}^2\text{d}x}+\frac{\rho }{2}{\displaystyle {\int }_{\Omega }{\left|I_0-B-S^k\right|}^2\text{d}x}+\langle {\theta }^k,I_0-B-S^k\rangle ,\end{array}$ (15)

对应的欧拉–拉格朗日方程为

$\rho B^{k+1}-2\alpha \Delta B^{k+1}={\theta }^k+\rho \left(I_0-S^k\right),$ (16)

利用快速傅里叶变换，可得

$B^{k+1}={ℱ}^{-1}\left(\frac{ℱ\left[\rho \left(I_0-S^k\right)+{\theta }^k\right]}{\rho -2\alpha ℱ\left(\Delta \right)}\right).$ (17)

② 固定$u^k,S^k,d_1^{k+1},d_2^{k+1},B^{k+1},{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k$ ，求解$W^{k+1}$ 。

$\begin{array}{c}{W}^{k+1}=\underset{W}{\arg \min }L\left(u^k,S^k,B^{k+1},{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k,d_1^{k+1},d_2^{k+1}\right)\\ =\underset{W}{\arg \min }{\lambda }_1{\displaystyle {\int }_{\Omega }{\left(c_1-S^k\right)}^{2}W\text{d}x}+\frac{{\rho }_3}{2}{\displaystyle {\int }_{\Omega }{\left|W-u^k\right|}^2}\text{d}x\\ +{\lambda }_2{\displaystyle {\int }_{\Omega }{\left(c_2-S^k\right)}^2}\left(1-W\right)\text{d}x+\langle {\theta }_3,W-u^k\rangle ,\end{array}$ (18)

求导得

$W^{k+1}=\mathrm{Pr}oj\frac{{\lambda }_2{\left(c_2-S^k\right)}^2-{\lambda }_1{\left(c_1-S^k\right)}^2-{\theta }_3{}^k+{\rho }_3{u}^k}{{\rho }_3},$ (19)

其中$pro{j}_{\Gamma }\left(W^{k+1}\right):=\max \left\{\min \left\{W^{k+1},1\right\},0\right\}$ 。

③ 固定$u^k,W^{k+1},d_1^{k+1},d_2^{k+1},B^{k+1},{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k$ ，求解$S^{k+1}$ 。

$\begin{array}{c}{S}^{k+1}=\underset{S}{\arg \min }L\left(u^k,W^{k+1},B^{k+1},{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k,d_1^{k+1},d_2^{k+1}\right)\\ =\underset{S}{\arg \min }{\lambda }_1{\displaystyle {\int }_{\Omega }{\left(c_1-S\right)}^2{W}^{k+1}\text{d}x}+{\lambda }_2{\displaystyle {\int }_{\Omega }{\left(c_2-S\right)}^2}\left(1-W^{k+1}\right)\text{d}x\\ +\frac{\rho }{2}{\displaystyle {\int }_{\Omega }{\left|I_0-B^{k+1}-S\right|}^2}\text{d}x+\langle {\theta }^k,I_0-B^{k+1}-S\rangle \\ +\frac{{\rho }_2}{2}{\displaystyle {\int }_{\Omega }{\left|d_2^{k+1}-\nabla S\right|}^2}\text{d}x+\langle {\theta }_2^k,\nabla S-d_2^{k+1}\rangle ,\end{array}$ (20)

对应的欧拉–拉格朗日方程为

$\begin{array}{l}\left(\rho +2{\lambda }_1{u}^{k+1}+2{\lambda }_2-2{\lambda }_2{u}^{k+1}\right)S^{k+1}-{\rho }_2\Delta S^{k+1}\\ =\nabla \cdot \left({\theta }_2^k-d_2^{k+1}{\rho }_2\right)+\rho \left(I_0-B^{k+1}\right)+{\theta }^{k+1}+2c_1{\lambda }_1{u}^{k+1}+2c_2{\lambda }_2-2c_2{\lambda }_2{u}^{k+1},\end{array}$ (21)

在后续的数值实验中，保持参数${\lambda }_1$ 和${\lambda }_2$ 相等。因此这里用快速傅里叶变换进一步更新$S^{k+1}$ ，

$S^{k+1}={ℱ}^{-1}\left(\frac{ℱ\left[\nabla \cdot \left({\theta }_2^k-d_2^{k+1}{\rho }_2\right)+\rho \left(I_0-B^{k+1}\right)+{\theta }^{k+1}+f^{k+1}\right]}{\rho +2{\lambda }_2-{\rho }_2ℱ\left(\Delta \right)}\right),$ (22)

其中$f^{k+1}=2c_1{\lambda }_1{u}^{k+1}+2c_2{\lambda }_2-2c_2{\lambda }_2{u}^{k+1}$ 。

④ 固定$S^{k+1},W^{k+1},d_1^{k+1},d_2^{k+1},B^{k+1},{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k$ ，求解$u^{k+1}$ 。

$\begin{array}{c}{u}^{k+1}=\underset{u}{\arg \min }L\left(W^{k+1},S^{k+1},B^{k+1},{\theta }^k,{\theta }_1^k,{\theta }_2^k,{\theta }_3^k,d_1^{k+1},d_2^{k+1}\right)\\ =\underset{u}{\arg \min }\frac{{\lambda }_3}{2}u\left(u-f\left(u\right)\right)+\frac{{\rho }_1}{2}{\displaystyle {\int }_{\Omega }{\left|d_1^{k+1}-\nabla u\right|}^2}\text{d}x+\langle {\theta }_1^k,\nabla u-d_1^{k+1}\rangle \\ +\frac{{\rho }_3}{2}{\displaystyle {\int }_{\Omega }{\left|W^{k+1}-u\right|}^2}\text{d}x+\langle {\theta }_3^k,W^{k+1}-u\rangle ,\end{array}$ (23)

等价于求解

$\underset{u}{\arg \min }\frac{{\rho }_1}{2}{\Vert \nabla u-d_1^{k+1}+\frac{{\theta }_1^k}{{\rho }_1}\Vert }_2^2+\frac{{\rho }_3}{2}{\Vert W^{k+1}-u+\frac{{\theta }_3^k}{{\rho }_3}\Vert }_2^2+\frac{{\lambda }_3}{2}u\left(u-f\left(u\right)\right),$ (24)

根据欧拉–拉格朗日方程，可以得到

${\rho }_{3}u-{\rho }_1△u+{\lambda }_{3}u={\lambda }_{3}f\left(u\right)+{\rho }_3{W}^{k+1}-{\theta }_3^k-{\rho }_{1}div{d}_1^{k+1}+div{\theta }_1^k,$ (25)

利用快速傅里叶变换，

$u^{k+1}={ℱ}^{-1}\left\{\frac{ℱ\left({\lambda }_{3}f\left(u\right)+{\rho }_3{W}^{k+1}-{\theta }_3^k-{\rho }_{1}div{d}_1^{k+1}+div{\theta }_1^k\right)}{{\rho }_3-{\rho }_1ℱ\left(\Delta \right)+{\lambda }_3}\right\}.$ (26)

其中本文所使用的去噪器$f\left(u\right)$ 为FFDNet [23]。

⑤ 固定$u^{k+1},S^{k+1},W^{k+1},B^{k+1},d_1^{k+1},d_2^{k+1}$ ，求解${\theta }^{k+1},{\theta }_1^{k+1},{\theta }_2^{k+1},{\theta }_3^{k+1}$ 。

${\theta }^{k+1}={\theta }^k+\rho \left(I_0-B^{k+1}-S^{k+1}\right),$ (27)

${\theta }_1^{k+1}={\theta }_1^k+{\rho }_1\left(\nabla u^{k+1}-d_1^{k+1}\right),$ (28)

${\theta }_2^{k+1}={\theta }_2^k+{\rho }_2\left(\nabla S^{k+1}-d_2^{k+1}\right),$ (29)

${\theta }_3^{k+1}={\theta }_3^k+{\rho }_3\left(W^{k+1}-u^{k+1}\right).$ (30)

综合前面的结果，变分问题(8)的具体算法步骤如下。

步骤1 $S^0=I^0,B^0=\left(0,0\right),W^0=0,d_1^0=d_2^0=\left(0,0\right),{\theta }^0={\theta }_1^0={\theta }_2^0={\theta }_3^0\left(0,0\right)$ 。

步骤2 迭代更新k步。

$d_1^{k+1}=\frac{\nabla u^k+\frac{{\theta }_1^k}{{\rho }_1}}{\left|\nabla u^k+\frac{{\theta }_1^k}{{\rho }_1}\right|}\max \left\{\left|\nabla u^k+\frac{{\theta }_1^k}{{\rho }_1}\right|-\frac{1}{{\rho }_1},0\right\},$

$d_2^{k+1}=\frac{\nabla S^k+\frac{{\theta }_2^k}{{\rho }_2}}{\left|\nabla S^k+\frac{{\theta }_2^k}{{\rho }_2}\right|}\max \left\{\left|\nabla S^k+\frac{{\theta }_2^k}{{\rho }_2}\right|-\frac{\beta }{{\rho }_2},0\right\},$

$B^{k+1}={ℱ}^{-1}\left(\frac{ℱ\left[\rho \left(I_0-S^k\right)+{\theta }^k\right]}{\rho -2\alpha ℱ\left(\Delta \right)}\right),$

$W^{k+1}=Proj\frac{{\lambda }_2{\left(c_2-S^k\right)}^2-{\lambda }_1{\left(c_1-S^k\right)}^2-{\theta }_3^k+{\rho }_3{u}^k}{{\rho }_3},$

$S^{k+1}={ℱ}^{-1}\left(\frac{ℱ\left[\nabla \cdot \left({\theta }_2^k-d_2^{k+1}{\rho }_2\right)+\rho \left(I_0-B^{k+1}\right)+{\theta }^{k+1}+f^{k+1}\right]}{\rho +2{\lambda }_2-{\rho }_2ℱ\left(\Delta \right)}\right),$

$u^{k+1}={ℱ}^{-1}\left\{\frac{ℱ\left({\lambda }_{3}f\left(u\right)+{\rho }_3{W}^{k+1}-{\theta }_3^k-{\rho }_{1}div{d}_1^{k+1}+div{\theta }_1^k\right)}{{\rho }_3-{\rho }_1ℱ\left(\Delta \right)+{\lambda }_3}\right\},$

${\theta }^{k+1}={\theta }^k+\rho \left(I_0-B^{k+1}-S^{k+1}\right),$

${\theta }_1^{k+1}={\theta }_1^k+{\rho }_1\left(\nabla u^{k+1}-d_1^{k+1}\right),$

${\theta }_2^{k+1}={\theta }_2^k+{\rho }_2\left(\nabla S^{k+1}-d_2^{k+1}\right),$

${\theta }_3^{k+1}={\theta }_3^k+{\rho }_3\left(W^{k+1}-u^{k+1}\right).$

步骤3 设置迭代终止条件。

4. 实验结果及分析

本节展示数值实验，旨在说明本文模型和算法的有效性。主要展示纹理、结构和MR图像的分割结果，并与两种基于深度学习的分割方法：GateNet [26]、BASNet [27]以及三种传统模型的分割方法：CVBS模型[11]、LCVB模型[12]和AITV_PS模型[28]进行了比较。

使用Dice相似系数(DSC) [29]和Hausdorff距离(HD) [30]作为不同算法分割准确性的定量评估指标。HD是表示最显著分割误差的指标，DSC是衡量分割结果与标准分割结果之间相似度度量，较小的HD和较大的DSC表示更精确地分割结果。

4.1. 参数设置及相关说明

本文模型设置参数${\lambda }_1={\lambda }_2$ ，$\beta =3$ ，${\lambda }_3$ 与输入FFDNet [25]噪声级图有关，它是估计值与标准分割结果差级相关之间的误差，然后根据偏置场调整平滑参数$\alpha$ 。具体参数设置如表1所示。在算法中，观察到所提出算法的数值性能对$\rho ,{\rho }_1,{\rho }_2$ 和${\rho }_3$ 的敏感性相对较低，只要这些值不落入过大或过小的极端值。

因此，本文设置$\rho =50,{\rho }_1=0.1,{\rho }_2=0.1$ 和${\rho }_3=0.2$ ，并且建立了停止准则：$\frac{\Vert u^{k+1}-u^k\Vert }{\Vert u^k\Vert }\le ϵ$ 。

Table 1. Parameters of all pictures in the experiment

表1. 实验中所有图片的参数

图像	${\lambda }_k$	$\alpha$	图像	${\lambda }_k$	$\alpha$
1	${\lambda }_1={\lambda }_2=105,{\lambda }_3=0.1$	20	9	${\lambda }_1={\lambda }_2=55,{\lambda }_3=0.2$	50
2	${\lambda }_1={\lambda }_2=30,{\lambda }_3=0.2$	100	10	${\lambda }_1={\lambda }_2=25,{\lambda }_3=0.2$	10
3	${\lambda }_1={\lambda }_2=5,{\lambda }_3=0.2$	500	11	${\lambda }_1={\lambda }_2=80,{\lambda }_3=0.1$	100
4	${\lambda }_1={\lambda }_2=105,{\lambda }_3=0.2$	500	12	${\lambda }_1={\lambda }_2=100,{\lambda }_3=0.1$	100
5	${\lambda }_1={\lambda }_2=10,{\lambda }_3=0.2$	5000	13	${\lambda }_1={\lambda }_2=4,{\lambda }_3=0.1$	100
6	${\lambda }_1={\lambda }_2=30,{\lambda }_3=0.2$	100	14	${\lambda }_1={\lambda }_2=5,{\lambda }_3=0.1$	100
7	${\lambda }_1={\lambda }_2=20,{\lambda }_3=0.2$	500	15	${\lambda }_1={\lambda }_2=7,{\lambda }_3=0.1$	100
8	${\lambda }_1={\lambda }_2=5,{\lambda }_3=0.2$	500	16	${\lambda }_1={\lambda }_2=6,{\lambda }_3=0.1$	500

为了验证模型的分割能力，将所提出的模型与其他对比模型在Weizmann分割数据集以及脑磁共振数据集[31]的所有图像上进行比较，发现本文模型在这两个数据集中几乎所有图片都能表现出最好的效果。为了方便说明，本文选取的10幅自然图像都能明显地看出丰富的纹理特征以及轮廓鲜明的结构，4张MR图像灰质和白质线条分明，对于分割结果可以很容易地观察出效果。

4.2. 初始轮廓鲁棒性实验

本小节检验提出的方法对初始轮廓的鲁棒性。图2测试了具有不同初始轮廓的简单合成图像的性能。图2(a)代表输入图像，图2(b)对应地面真值。图2(c)~(f)描绘了四个不同的初始轮廓。图2(c)的初始化包含整个五角星，图2(d)和图2(e)的初始化分为五角星两侧的交界处和五角星的一个顶点，图2(f)的初始化完全在五角星内部。这些初始化方法差别很大，但四种方法分割结果的HD值都是0，这说明提出的方法对轮廓初始化具有较好地鲁棒性。

(a) (b)

(c) (d) (e) (f)

(g) (h) (i) (j)

Figure 2. Segmentation results of different initial contours. (a) Original image; (b) Ground truth; (c) Initial outline 1; (d) Initial outline 2; (e) Initial outline 3; (f) Initial outline 4; (g) Result of contour 1; (h) Result of contour 2; (i) Result of contour 3; (j) Result of contour 4

图2. 不同初始轮廓的分割结果。(a) 原始图像；(b) 标准分割结果；(c)初始轮廓1；(d) 初始轮廓2；(e) 初始轮廓3；(f) 初始轮廓4；(g) 轮廓1分割结果；(h) 轮廓2分割结果；(i) 轮廓3分割结果；(j) 轮廓4分割结果

4.3. 纹理图像分割

本节对六幅纹理图像(见图3)进行了检验，以证明所提出方法的有效性。实验结果如图4所示，从图中可以看出，GateNet模型和BASNet模型在对图像1~6分割结果都与标准分割结果相差很大，不能反映图片的实际情况；CVBS模型和LCVB模型在处理纹理图片时，二者的分割效果差不多，但是从图像1可以看出，这两种方法出现了比较严重的过分割问题，花蕊实际上不需要划分出来，但是它们的分割结果却保留了大部分的花蕊部分。从图像4中树叶内部又可以明显地看出这两种方法存在着欠分割的问题；AITV_PS模型总体上对于纹理图像的分割效果较好，但是图像2铁塔的底部、图像3铜锁的标签都还有一定的瑕疵。与这些方法相比，提出的模型有效地融合了耦合TV捕获边缘和CNN捕获细节的优势，减少了欠分割、过分割等问题，对图像的细节和纹理部分具有更好的分割效果。

(a) (b) (c) (d) (e) (f)

Figure 3. Original texture images. (a) Image 1; (b) Image 2; (c) Image 3; (d) Image 4; (e) Image 5; (f) Image 6

图3. 原始纹理图像。(a) 图像1；(b) 图像2；(c) 图像3；(d) 图像4；(e) 图像5；(f) 图像6

(a) (b) (c) (d) (e) (f) (g)

Figure 4. Results of texture image segmentation using different models and standard segmentation results. (a) CVBS model; (b) LCVB model; (c) GateNet model; (d) BASNet model; (e) AITV_PS model; (f) The proposed model; (g) Ground truth

图4. 使用不同模型得到的纹理图像分割结果及标准分割结果。(a) CVBS模型；(b) LCVB模型；(c) GateNet模型；(d) BASNet模型；(e) AITV_PS模型；(f) 本文模型；(g) 标准分割结果

为了进一步验证纹理图像的分割能力，将本文方法的分割结果与其他模型的分割结果进行定量比较。各模型对应的HD和DSC值如表2所示。实验结果表明，该方法对细节丰富的纹理图像的分割效果最好。

Table 2. DSC values and HD values of the segmentation results in Figure 4

表2. 图4中分割结果的DSC值和HD值

模型		CVBS	LCVB	GateNet	BASNet	AITV_PS	Ours
1	DSC	0.9762	0.9785	0.9632	0.9692	0.9819	0.9939
	HD	1.4114	1.3268	1.8019	1.6486	1.2012	0.6479
2	DSC	0.9733	0.9793	0.9536	0.9663	0.9791	0.9872
	HD	1.5661	1.2710	2.2202	1.8296	1.1707	0.9275
3	DSC	0.9843	0.9831	0.9536	0.9588	0.9857	0.9908
	HD	0.5694	0.5328	0.8781	0.8369	0.5169	0.3871
4	DSC	0.9626	0.9686	0.9425	0.9630	0.9738	0.9824
	HD	2.0877	1.8866	2.8088	2.0609	1.7270	1.3226
5	DSC	0.9715	0.9747	0.9586	0.9678	0.9729	0.9785
	HD	0.9604	0.8800	1.1034	0.9841	0.9069	0.8546
6	DSC	0.9741	0.9746	0.9521	0.9686	0.9677	0.9777
	HD	0.6812	0.6418	0.9548	0.7287	0.7477	0.6194

4.4. 结构图像分割

本节展示每种分割方法对结构图像(见图5)的分割效果。实验结果如图6所示。从图6可以看出，两种深度学习的方法分割出的结果与标准结果相差很大，这是由于深度学习方法是通过端到端的方式从原始数据中学习特征表示，如果特征学习不充分，会导致分割结果不佳；从图像7和11可以看出，CVBS模型出现了边界模糊的问题；从图像10可以看出，三种传统方法都有一定的边界丢失；AITV_PS模型对于结构图像分割结构效果较好，但是对于图像9中十字架的尖角特征，图像12小狗的腿部轮廓并没有很好地体现出来。从本文提出的方法对这六幅图像的分割结果来看，本文的方法使用深度图像先验并在灵活的RED框架下捕获了更多的细节，因此分割出来的物体更准确。

(a) (b) (c) (d) (e) (f)

Figure 5. Original structure images. (a) Image 7; (b) Image 8; (c) Image 9; (d) Image 10; (e) Image 11; (f) Image 12

图5. 原始结构图像。(a) 图像7；(b) 图像8；(c) 图像9；(d) 图像10；(e) 图像11；(f) 图像12

(a) (b) (c) (d) (e) (f) (g)

Figure 6. Results of structure image segmentation using different models and standard segmentation results. (a) CVBS model; (b) LCVB model; (c) GateNet model; (d) BASNet model; (e) AITV_PS model; (f) The proposed model; (g) Ground truth

图6. 使用不同模型得到的结构图像分割结果及标准分割结果。(a) CVBS模型；(b) LCVB模型；(c) GateNet模型；(d) BASNet模型；(e) AITV_PS模型；(f) 本文模型；(g) 标准分割结果

为了客观评价分割结果的准确性，表3计算并给出了它们对应的HD和DSC值。提出的方法具有最高的DSC值和最低的HD值，这证明了本文模型对结构图像的分割效果最好。

Table 3. DSC values and HD values of the segmentation results in Figure 6

表3. 图6中分割结果的DSC值和HD值

模型		CVBS	LCVB	GateNet	BASNet	AITV_PS	Ours
7	DSC	0.9688	0.9649	0.9794	0.9461	0.9691	0.9929
	HD	1.7859	1.5721	0.9666	2.0308	1.4860	0.6671
8	DSC	0.9734	0.9757	0.9647	0.9681	0.9819	0.9889
	HD	1.5652	1.2624	1.7033	1.5408	1.0995	1.0871
9	DSC	0.9582	0.9533	0.9032	0.9101	0.9612	0.9696
	HD	1.0378	0.9095	1.6458	1.5860	0.8942	0.8427
10	DSC	0.9770	0.9783	0.9539	0.9689	0.9831	0.9883
	HD	0.8921	0.8888	1.1497	0.9927	0.8235	0.7001
11	DSC	0.8057	0.8388	0.9043	0.8666	0.8863	0.9406
	HD	3.8300	2.8123	2.4049	3.2373	2.8268	1.6072
12	DSC	0.8125	0.9152	0.9238	0.9210	0.9215	0.9728
	HD	3.2210	2.0321	1.7661	1.7792	1.7693	0.1224

4.5. MR图像分割

本节介绍了各种模型在脑磁共振图像(如图7)分割的性能，分割结果如图8所示。BASNet和GateNet两种深度学习模型对MR图像的分割效果十分不理想，比起上两小节对结构图像和纹理图像处理的结果情况更糟糕，其主要原因是这两种网络是针对自然图像分割设计的，而对于医学图像会更容易地出现特征学习不足的情况；三种传统模型对于MR图像的分割结果都比较精确，但是在图中绿框区域，提出的模型可以更好的捕捉到这些区域的细节，使得分割结果比这三种模型更加细致。

(a) (b) (c) (d)

Figure 7. MR images. (a) Image 13; (b) Image 14; (c) Image 15; (d) Image 16

图7. MR图像。(a) 图像13；(b) 图像14；(c) 图像15；(d) 图像16

Figure 8. Results of MR image segmentation using different models and standard segmentation results. (a) CVBS model; (b) LCVB model; (c) GateNet model; (d) BASNet model; (e) AITV_PS model; (f) The proposed model; (g) Ground truth

图8. 使用不同模型得到的MR图像分割结果及标准分割结果。(a) CVBS模型；(b) LCVB模型；(c) GateNet模型；(d) BASNet模型；(e) AITV_PS模型；(f) 本文模型；(g) 标准分割结果

此外，表4给出了各种方法得到的分割结果所对应的HD和DSC值。从表中可以明显看出，提出的方法拥有最低的HD值和最高的DSC值。这些结果证实了提出的方法在MR脑图像分割中具有较高的准确性。

Table 4. DSC values and HD values of the segmentation results in Figure 8

表4. 图8中分割结果的DSC值和HD值

模型		CVBS	LCVB	GateNet	BASNet	AITV_PS	Ours
13	DSC	0.8884	0.8880	0.8485	0.8651	0.8992	0.9024
	HD	2.8588	2.8815	3.3677	3.1540	2.7465	2.6455
14	DSC	0.8725	0.9039	0.8287	0.8448	0.9037	0.9101
	HD	3.1050	2.6402	3.5619	3.4012	2.4428	2.5996
15	DSC	0.7911	0.7967	0.7080	0.7341	0.7991	0.8193
	HD	2.6706	2.7229	3.3580	3.1985	2.6999	2.5490
16	DSC	0.8921	0.8987	0.8570	0.8710	0.9095	0.9148
	HD	2.9169	2.8235	3.3860	3.1969	2.6800	2.5154

4.6. 计算时间

在相同设备上，本文将两个数据集中共120张图像分别在CVBS模型[11]、LCVB模型[12]、GateNet模型[26]、BASNet模型[27]、AITV_PS模型[28]以及本文模型上进行运算，各个模型的运算时间如表5所示。从表5可以看出，两种深度学习的方法在运行时间上最快，这是因为深度学习模型通常使用图形处理单元(Graphics Processing Unit, GPU)进行计算，而GPU具有强大的并行计算能力，能够同时处理大量数据和计算任务，从而加快模型的运行速度。本文模型在CVBS模型上加入去噪正则项，虽然模型复杂度变高，但是由于RED灵活高效以及融入深度先验的优势，使其在运行速度上并没有落后CVBS模型很多。此外，对于LCVB模型和AITV_PS模型，本文模型的运行效率显著高于二者。综上所述，本文模型具有快速高效的优势。

Table 5. Average operation time of images

表5. 图像的平均运算时间

模型	CVBS	LCVB	GateNet	BASNet	AITV_PS	Ours
平均运行时间	230.25	750.36	40.21	63.85	195.66	105.24

5. 结语

本文从CVBS模型[11]出发，结合CNN捕捉细节的特征，提出了一种融合深度图像先验的TV分割模型。在去噪正则化框架下，采用ADMM算法对模型进行数值求解，并证明了模型的鲁棒性。通过与CVBS模型[11]、LCVB模型[12]、GateNet模型[26]、BASNet模型[27]和AITV_PS模型[28]的分割结果比较，发现该模型有效地耦合了TV捕获边缘和CNN捕获细节的功能，对结构纹理和MR图像的分割效果显著。此外，计算出各个模型的运行时间，分析得出本文模型具有算法复杂度低，运行速度快的优势。综上所述，本文模型在处理结构丰富和纹理细致的图像时，可以分割出更多的细节，分割结果更精确，并且结合深度神经网络，在运行效率上相比于传统模型显著提升。

NOTES

^*通讯作者。