发布时间: 2020-04-16
摘要点击次数:
全文下载次数:
DOI: 10.11834/jig.180659
2020 | Volume 25 | Number 4

医学图像处理

低秩和主动形状模型的肺区分割

孙申申¹, 田丹¹, 吴微¹, 康雁², 赵宏²

1. 沈阳大学信息工程学院, 沈阳 110044;

2. 东北大学中荷生物医学与信息工程学院, 沈阳 110016

收稿日期: 2019-01-14; 修回日期: 2019-09-20; 预印本日期: 2019-09-27

基金项目: 国家自然科学基金项目(61703285);中国博士后科学基金面上项目(2019M651142);辽宁省自然科学基金指导计划项目(20180550615);沈阳市科技计划项目(18013015)

第一作者简介: 孙申申, 1980年生, 女, 副教授, 硕士生导师, 主要研究方向为计算机辅助诊断。E-mail:sunshenshen@126.com;
田丹, 女, 副教授, 硕士生导师, 主要研究方向为图像处理。E-mail:www.sltd2008@163.com;
吴微, 女, 副教授, 硕士生导师, 主要研究方向为图像处理。E-mail:wuwei429@163.com;
康雁, 男, 教授, 主要研究方向为计算机辅助诊断。E-mail:Kangyan@SZtu.edu.cn;
赵宏, 男, 教授, 主要研究方向为计算机辅助诊断。E-mail:zhaoh@neusoft.com.

中图法分类号: TP391.4

文献标识码: A

文章编号: 1006-8961(2020)04-0759-09

摘要

目的肺区分割是肺癌计算机辅助诊断系统的首要步骤。主动形状模型（active shape model，ASM）能根据训练集获得肺区形状模型，再结合待分割肺区影像自身的局部特征，进行测试影像的分割。由于主成分分析（principal component analysis，PCA）仅能去除服从高斯分布的噪声，不能处理其他类型的噪声，所以当训练集含有非高斯类型的噪声样本时，采用基于PCA的ASM无法训练出正确的形状模型，使得肺区分割不能得到正确的结果。而低秩（low rank，LR）理论的鲁棒主成分分析（robust principal component analysis，RPCA）能去除各种类型的噪声，基于此，本文提出一种将RPCA与ASM相结合的方法。方法首先对训练样本集标记点矩阵进行低秩分解，去除噪声样本对训练出的形状模型的影响。然后在ASM训练局部梯度模型时，用判断训练样本轮廓上的标记点曲率直方图的相似度来去除噪声样本。结果在训练集含噪声样本时，将基于RPCA的ASM与传统ASM（即基于PCA的ASM）分别生成的形状模型进行对比，发现基于RPCA的ASM生成的形状模型与训练集无噪声样本时传统ASM生成的形状模型更相符。在训练集含噪声样本的情况下，基于RPCA的ASM方法分割EMPIRE10数据集中的22个肺影像，与金标准的重叠度为94.5%，而基于PCA的ASM方法分割准确率仅为69.5%。结论实验结果表明，在训练样本集中有噪声样本的情况下，基于RPCA的ASM分割能得到比基于PCA的ASM更好的分割效果。

关键词

低秩; 主动形状模型; 鲁棒主成分分析; 肺区分割; 噪声样本

Lung segmentation by active shape model approach based on low rank theory

Sun Shenshen¹, Tian Dan¹, Wu Wei¹, Kang Yan², Zhao Hong²

1. Information and Engineering School, Shenyang University, Shenyang 110044, China;

2. Sino-Dutch Biomedical and Information Engineering School, Northeastern University, Shenyang 110016, China

Supported by: National Natural Science Foundation of China (61703285)

Abstract

Objective Computer-aided diagnosis of lung cancer is the automatic recognition of lung lesions in computed tomography (CT) images through a computer image processing technology. Given that lung lesions are located in the lung areas, the contours of lung areas should be automatically marked on lung images first. Thus, lung segmentation is the first step in the computer-aided diagnosis of lung cancer. Current segmentation methods for pulmonary parenchyma, including normal and pathological, can be divided into three categories according to visual appearance, location relationship of shape, and whether or not the method is a hybrid. Methods based on visual appearance mostly consider the lower gray values of lung areas instead of surrounding tissues. Therefore, a pathological lung (especially including interstitial lung disease) is segmented by a classifier on the basis of gray and texture features. These methods mainly introduce the anatomical relationship of the lungs with the heart, liver, spleen, and rib. Hybrid methods establishes shape models and evolves lung contours according to surface features. They can achieve good segmentation results without noise samples in training sets but does not consider the problem of noise samples in the training sets during the construction of training shape models. The active shape model (ASM) based on principal component analysis (PCA) is an outstanding method, which can only deal with Gaussian noise. ASM is a statistical deformation model and uses PCA to obtain the average shape and the allowable range of shape of a training sample set and then establishes a shape model. ASM establishes a local feature appearance model according to the training sample set. Finally, ASM evolves the contour of the test image on the basis of the local feature appearance and shape models. Robust PCA (RPCA) in low rank theory is robust to any kind of noise. In addition, the shapes of a training sample set have a low rank attribute because the anatomical shapes of the human lungs are roughly the same. Therefore, this study combined RPCA with ASM to solve the problem that the correct shape model cannot be obtained because of noise samples in training sets. Method First, a traditional method was used to segment the pulmonary region of training samples and generate lung boundary contours. The segmentation steps include extra-pulmonary region exclusion by iteration threshold and region growing segmentation methods, trachea elimination, internal cavity filling, liver elimination, and left and right lung separation. Second, we extracted 96 marker points on the left or right lung boundary contour. We assigned these marker points as $x $-coordinates and $y $-coordinates and placed them into a single matrix. Then, we aligned the lung boundary contours in the training sample set through the Procrustes analysis method. Third, we decomposed the matrix through the low rank theory into the low rank component and noise component matrices. We used the accelerated approximation gradient method to solve the decomposition and employed PCA to the low rank component matrix and obtained the eigenvalues and eigenvectors. They were used to construct the shape model, which is expressed as the average shape and deformation range. Then, we established a local feature appearance model by using gradient features and obtained a gradient-level appearance model that describes the typical image structure around each landmark from pixel profiles sampled (using linear interpolation) around each landmark perpendicular to the contour. Finally, ASM evolved the contour of the test image based on the local feature appearance and shape models. In the training of the local gradient model in ASM, the noise samples of the training set must be eliminated. In this study, the similarity of the curvature histograms of the marked points on the contours of the training samples was evaluated, and noise samples were removed. Result First, a training dataset was established. It consisted of 35 multi-slice spiral CT lung scan images collected from the First Affiliated Hospital of Guangzhou Medical University. Two scans included juxta-pleural tumor on the left side. The test dataset comprised 22 scans from EMPIRE 10 and is the gold standard for lung contours.The five scans included juxta-pleural tumor on the left side. Using the proposed method to segment the 22 lungs of EMPIRE10 database, we obtained the overlap rate with a gold standard of 94.5%. Meanwhile, the accuracy of ASM based on PCA was only 69.5%. The shape models generated by RPCA-based ASM and traditional ASM (i.e., PCA-based ASM) were compared when the training set contained noise samples. Notably, the shape model generated by ASM based on RPCA was consistent with the shape model generated by the traditional ASM when the training set was noiseless. Remarkably, RPCA eliminated the influence of the noise samples on the shape model. Adding PRCA to other lung parenchyma segmentation algorithms generated good results for training sets containing noise samples. The gray threshold method and ASM segmentation method based on RPCA were used to segment the same lung image with juxta-pleural tumor.Notably, the ASM segmentation method based on RPCA generated correct segmentation results in contrast to the gray threshold method. Conclusion The experimental results showed that a good segmentation result was obtained by the ASM based on RPCA when the noise sample was contained in the training sample set. When the pulmonary parenchyma was absent, some markers were updated to the transition region of the normal lung tissues and tumor rather than to the real lung border during searching.Thus, identifying abnormal markers first is necessary. A new abnormal marker identification method should be designed to provide different searching functions for normal and abnormal markers.

Key words

low rank(LR); active shape model (ASM); robust principal component analysis (RPCA); lung segmentation; noise samples

0 引言

肺癌计算机辅助诊断是利用计算机图像处理技术对肺部CT影像中的病灶进行自动识别。由于病灶长在肺区内，所以首先应在肺影像上自动标出肺区的轮廓，故肺区分割是肺癌计算机辅助诊断的首要步骤。

正常肺区和病态肺区(即肺内含有大面积的病灶)分割的方法大致分为3类：1)基于表面灰度或纹理信息的分割方法(Hosseini-Asl等，2016)，大多认为肺区的灰度值比周围其他组织低，故用基于灰度、纹理等特征的分类器进行病态肺(特别是含有间质性肺疾病)的肺区分割；2)基于解剖组织位置关系的分割方法(Oluyide等，2018)，考虑肺区与心、肝脏、脾和肋骨之间的解剖位置关系；3)基于形状模型(Soliman等，2017)的分割方法，用肺区和纵隔区的局部和全局表面特征建立形状模型。例如，Sun等人(2012)通过主动形状模型(active shape model，ASM)得到肺初始边界，然后用全局表面特征驱动轮廓的演化；Sluimer等人(2005)用15个肺影像序列建立正常肺区概率图，并用配准的方法做病态肺分割；Sofka等人(2011)用在肋骨和脊柱上的标记点初始化形状模型，再结合可变形模型分割含有大肿块的肺区。但上述方法在训练形状模型的过程中都没有考虑训练集中含有噪声样本的问题。

ASM是一种基于统计的变形模型(Cootes等，1995)，使用主成分分析(principal component analysis，PCA)获得训练样本集的平均形状和形状允许的变化范围，从而建立形状模型。目标图像结合自身特征在形状模型的驱动下进行分割(van Ginneken等，2002)。但ASM中使用的PCA仅能去除服从高斯分布的噪声样本的影响，不能处理其他类型的噪声样本。即当训练集里存在非高斯类型的噪声样本时，ASM无法训练出正确的形状模型(Zhou等，2014)。

由于训练样本集合的形状往往是低秩(low rank，LR)的，且低秩理论中鲁棒主成分分析(robust principal component analysis，RPCA)对噪声有更高的容忍度(Zhou和Yu，2013)。故本文将RPCA与ASM相结合，解决在训练集中存在噪声样本的情况下，训练不出正确的形状模型的问题。

不少研究将低秩理论用于图像形状处理(包括对齐和分割)。Peng等人(2012)将低秩理论用于图像集合的非刚性对齐问题。Yang等人(2014)提出固定可变形脸集合对齐算法，将低秩理论和分段仿射变换结合到对齐模型中。Cheng(2013)将低秩理论与主动表面模型(active appearance models，AAM)相结合，对形状和纹理表面变化同时加以低秩的限制。Zhang等人(2012)和Wang等人(2013)提出基于稀疏形状组合模型(sparse shape composition，SSC)的肝脏形状建模。以上方法解决了统计形状模型中的非高斯误差敏感性问题，取得了较好效果，但目前还没有将低秩理论应用到ASM模型的工作。

本文提出将低秩理论的RPCA与ASM相结合，即首先对训练样本集标记点矩阵和协方差矩阵进行低秩分解，去除噪声样本对训练出的形状模型的影响。在ASM训练局部梯度模型中，本文用判断训练样本轮廓上的标记点曲率直方图的相似度来去除噪声样本的影响。

1 相关工作

1.1 主成分分析

对矩阵进行主成分分析的目的是找到最能表达该矩阵的特征向量。将图像表示成向量并排列起来组成矩阵$\boldsymbol{D} \in {\bf{R}}^{m \times n}, \boldsymbol{D}=\boldsymbol{X}+\boldsymbol{E}, \boldsymbol{X}$是低秩矩阵，$\boldsymbol{E}$是噪声矩阵，$\operatorname{rank}(\boldsymbol{D}) < \min (m, n)$，则

$ \begin{aligned} &\min\limits_{X}\|\boldsymbol{D}-\boldsymbol{X}\|_{\mathrm{F}}^{2}\\ &\text { s. } t \text {, } \operatorname{rank}(X) \leqslant k \end{aligned} $

(1)

式中，$\|\cdot\|_{\mathrm{F}}$为F范数，即$\|\boldsymbol{Y}\|_{\mathrm{F}}=\sqrt{\sum\limits_{i, j} y_{i j}^{2}}$。用式(1)低秩分解$\boldsymbol{D}$(邱枫等，2019)，基于奇异值分解(singular value decomposition，SVD)，式(1)的解为

$\boldsymbol{X}=\sum\limits_{i=1}^{k} \boldsymbol{\sigma}_{i} \boldsymbol{u}_{i} \boldsymbol{v}_{i}^{\mathrm{T}}$

(2)

式中，$\left\{\boldsymbol{u}_{i}\right\}, \left\{\boldsymbol{v}_{i}\right\}$和$\left\{\boldsymbol{\sigma}_{i}\right\}, i=1, \cdots, k$，分别为前$k$个左奇异向量、右奇异向量和$\boldsymbol{D}$的奇异值。

1.2 鲁棒主成分分析

由于PCA是用F范数表示的，所以PCA能容忍含有高斯噪声的数据，但对数据中的非高斯噪声敏感，对此，用低秩理论中的RPCA可以很好地解决。将矩阵$\boldsymbol{D}$分解成低秩成分$\boldsymbol{X}$和稀疏成分$\boldsymbol{E}$的和。求解时同时最小化$\boldsymbol{X}$的秩和$\boldsymbol{E}$的势，即

$ \begin{array}{c} \min\limits_{X, E}\|\boldsymbol{X}\|_{*}+\lambda\|\boldsymbol{E}\|_{1} \\ \text { s. t. } \boldsymbol{X}+\boldsymbol{E}=\boldsymbol{D} \end{array} $

(3)

式中，参数$λ$控制$\boldsymbol{X}$的秩，$λ$的选择应该依靠噪声的程度来确定。$\|\boldsymbol{X}\|_{*}$为核范数，$\|\boldsymbol{E}\|_1$为${\rm L}_1$范数。

1.3 ASM模型

ASM模型主要包括根据训练集建立的形状模型和表面局部特征模型，测试图像基于这两种模型搜索边界。

形状模型的建立步骤为：

1) 提取训练集中的$n$幅图像的边缘轮廓，记录在其上$m$个标记点的坐标且依次串联成向量

$ \boldsymbol{x}_{i}=\left[x_{i 0}, y_{i 0}, \cdots, x_{i j}, y_{i j}, \cdots, x_{i m}, y_{i m}\right]^{\mathrm{T}} $

式中，($x_{ij}$, $y_{ij}$)是第$i$幅图像的$j$个标记点的坐标，$i=1, 2, …, n，j=1, 2, …, m$。对齐形状向量后，形成形状向量集

$ \boldsymbol{S}=\left\{\boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \cdots, \boldsymbol{x}_{n}\right\} $

2) 取$\boldsymbol{S}$的前$t$个特征值对应的特征向量，记为$\boldsymbol{P}=\left[\boldsymbol{P}_{1}, \boldsymbol{P}_{2}, \cdots, \boldsymbol{P}_{t}\right]$。用$\boldsymbol{x}=\boldsymbol{\bar{x}}+\boldsymbol{P b}$建立形状模型，其中，$\boldsymbol{\bar{x}}$为均值向量，$\boldsymbol{P}$为特征向量矩阵，$\boldsymbol{b}$为形状参数向量。

3) 通过各标记点的灰度、梯度、纹理等局部特征的均值和协方差矩阵建立表面局部特征模型。以$\boldsymbol{\bar{x}}$为测试图像的初始轮廓，在形状模型允许的变化范围内，基于表面局部特征模型在测试图像上为各标记点搜索更合适的位置，迭代搜索过程，直至收敛(Wright等，2009)。

2 低秩的主动形状模型

2.1 低秩分解

ASM在建立形状模型时，如果训练集中含有噪声样本，就会产生错误的形状模型，影响分割结果。而噪声类型不仅有高斯噪声，还有其他类型的噪声，这就需要去除训练集中的噪声。而ASM中的PCA只能去除高斯噪声，无法去除其他类型的噪声。故本文将ASM中的PCA换为RPCA，即将训练集样本做低秩分解处理，同时将ASM步骤中的PCA得到的协方差矩阵也进行低秩分解处理。

用传统方法分割训练集肺图像(孙申申等，2016)，得到肺轮廓。设训练集肺轮廓为$\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n}, $，表示成向量形式为

$\boldsymbol{D}=\left(\boldsymbol{d}_{1}, \boldsymbol{d}_{2}, \cdots, \boldsymbol{d}_{n}\right)^{\mathrm{T}}$

式中，$\boldsymbol{d}_{i}=\left(x_{i j}, y_{i j}\right)$为标记点坐标。对$\boldsymbol{D}$做低秩处理，得

$ \begin{array}{l} \min\limits_{\boldsymbol{X, E}} \operatorname{rank}(\boldsymbol{X})+\lambda\|\boldsymbol{E}\|_{0} \\ \text { s. t. } \boldsymbol{X}+\boldsymbol{E}=\boldsymbol{D} \end{array} $

(4)

式中，$\boldsymbol{X}$为分解的低秩成分，$\boldsymbol{E}$为稀疏成分。式(4)是非凸函数，没有唯一解。所以用核范数代替低秩函数，用${\rm L}_1$范数代替${\rm L}_0$范数，其中${\rm L}_0$为非零元素的个数，使得式(4)为凸函数，即

$ \begin{aligned} &\min _{\boldsymbol{X, E}}\|\boldsymbol{X}\|_{*}+\lambda\|\boldsymbol{E}\|_{1}\\ &\text { s. t. }\|\boldsymbol{D}-\boldsymbol{X}-\boldsymbol{E}\|_{\mathrm{F}} \leqslant \delta \end{aligned} $

(5)

式中，$λ$为大于0的正则参数，$δ$为一个非常小的正数。用拉格朗日函数将式(5)写成式(6)的形式

$ L_{\mu}(\boldsymbol{X}, \boldsymbol{E})=\|\boldsymbol{X}\|_{*}+\lambda\|\boldsymbol{E}\|_{1}+\frac{\mu}{2}\|\boldsymbol{D}-\boldsymbol{X}-\boldsymbol{E}\|_{\mathrm{F}}^{2} $

(6)

式中，$μ$是一个小的正数，利用加速接近梯度算法(Wright等，2009)来求解。然后，求出均值向量为$\boldsymbol{m}_{\boldsymbol{x}}=\boldsymbol{E}\{\boldsymbol{X}\}$及协方差矩阵为

$ \boldsymbol{\varSigma}=\boldsymbol{E}\left\{\left(\boldsymbol{X}-\boldsymbol{m}_{x}\right)\left(\boldsymbol{X}-\boldsymbol{m}_{x}\right)^{\mathrm{T}}\right\} $

再对协方差矩阵做降噪的低秩分解，得

$ \begin{array}{c} \min\limits_{\boldsymbol{\varSigma, \gamma}} \operatorname{rank}(\boldsymbol{\varPhi})+\lambda\|\boldsymbol{\gamma}\|_{0} & \\ \text { s. t. } \boldsymbol{\varSigma} =\boldsymbol{\varPhi}+\boldsymbol{\gamma} \end{array} $

(7)

式中，$\boldsymbol{\varPhi}$为分解的低秩成分，$\boldsymbol{\gamma}$为稀疏成分。式(7)是非凸函数，所以用核范数代替rank，用${\rm L}_1$范数代替${\rm L}_0$范数，使得式(7)为凸函数，即

$ \begin{array}{c} \min\limits_{\boldsymbol{\varPhi}, \gamma}\|\boldsymbol{\varPhi}\|_{*}+\lambda\|\boldsymbol{\gamma}\|_{1}\\ \text { s.t. }\|\boldsymbol{\varSigma}-\boldsymbol{\varPhi}-\boldsymbol{\gamma}\|_{\mathrm{F}} \leqslant \delta \end{array} $

(8)

式(8)可以写为

$ L_{\mu}(\boldsymbol{\varPhi}, \boldsymbol{\gamma})=\|\boldsymbol{\varPhi}\|_{*}+\lambda\left\|_{\boldsymbol{\gamma}}\right\|_{1}+\frac{\mu}{2}\|\boldsymbol{\varSigma}-\boldsymbol{\varPhi}-\boldsymbol{\gamma}\|_{\mathrm{F}}^{2} $

(9)

用加速接近梯度算法(Wright等，2009)求解式(9)，用$\boldsymbol{X}=\boldsymbol{m}_{x}+\boldsymbol{\varPhi} \boldsymbol{b}$建立形状模型，$\boldsymbol{b}$为形状模型的形状参数向量。

最后，在建立局部特征模型时，为了避免受到训练样本集中噪声样本的影响，需删除噪声样本。噪声样本的识别方法如下：求训练样本集各样本$\boldsymbol{d}_{i}$和均值样本$\boldsymbol{m}_{x}$的曲率直方图，然后比较它们的巴氏距离，即

$ d\left(\boldsymbol{H}_{1}, \boldsymbol{H}_{2}\right)=\sqrt{1-\sum\limits_{i} \frac{\sqrt{\boldsymbol{H}_{1}(i) \times \boldsymbol{H}_{2}(i)}}{\sum\limits_{i} \boldsymbol{H}_{1}(i) \times \sum\limits_{i} \boldsymbol{H}_{2}(i)}} $

(10)

式中，$\boldsymbol{H}_1$, $\boldsymbol{H}_2$为直方图。当距离$d>Δ$时($Δ$为阈值)，说明是噪声样本。

以$\boldsymbol{m}_{x}$为测试图像的初始轮廓，在形状模型$\boldsymbol{X}$允许的变化范围内，基于表面局部特征模型在测试图像上为各标记点搜索轮廓边界。

2.2 模型求解

用加速接近梯度算法(Wright等，2009)求解式(5)和式(9)，其中，svd为奇异值分解函数。

1) 求解式(5)。输入$\boldsymbol{D} \in \bf{R}^{m \times n}$，权重$\lambda, \boldsymbol{X}_{0}, $$\boldsymbol{X}_{-1} \leftarrow 0, \boldsymbol{E}_{0}, \boldsymbol{E}_{-1} \leftarrow 0, t_{0}, t_{-1} \leftarrow 0, \mu_{0} \leftarrow 0.99 \times$$\|\boldsymbol{D}\|_{2, 2}, \bar{\mu} \leftarrow 10^{-5} \mu_{0}$。

当没有收敛时循环，否则退出循环。

$ \begin{array}{c} \tilde{\boldsymbol{X}}_{k} \leftarrow \boldsymbol{X}_{k}+\frac{t_{k-1}}{t_{k}}\left(\boldsymbol{X}_{k}-\boldsymbol{X}_{k-1}\right)\\ \widetilde{\boldsymbol{E}}_{k} \leftarrow \boldsymbol{E}_{k}+\frac{t_{k-1}}{t_{k}}\left(\boldsymbol{E}_{k}-\boldsymbol{E}_{k-1}\right)\\ \boldsymbol{Y}_{k}^{A} \leftarrow \tilde{\boldsymbol{X}}_{k}-\frac{1}{2}\left(\tilde{\boldsymbol{X}}_{k}+\tilde{\boldsymbol{E}}_{k}-D\right)\\ (\boldsymbol{U}, \boldsymbol{S}, \boldsymbol{V}) \leftarrow {svd}\left(\boldsymbol{Y}_{k}^{A}\right)\\ \boldsymbol{X}_{k+1} \leftarrow \boldsymbol{U}\left[\boldsymbol{S}-\frac{\mu}{2} \boldsymbol{I}\right]_{+} \boldsymbol{V}^{*}\\ \boldsymbol{Y}_{k}^{E} \leftarrow \widetilde{\boldsymbol{E}}_{k}-\frac{1}{2}\left(\tilde{\boldsymbol{X}}_{k}+\widetilde{\boldsymbol{E}}_{k}-\boldsymbol{D}\right)\\ \boldsymbol{E}_{k+1} \leftarrow \operatorname{sign}\left[\boldsymbol{Y}_{k}^{E}\right] \times\left[\left|\boldsymbol{Y}_{k}^{E}\right|-\frac{\lambda \mu}{2} 11\right]_{+} \\ t_{k+1} \leftarrow \frac{1+\sqrt{1+4 t_{k}^{2}}}{2}, \mu \leftarrow \max (0.9 \mu, \bar{\mu}) \end{array} $

循环结束，输出$\boldsymbol{X, E}$。

2) 求解式(9)。输入$\boldsymbol{\varSigma} \in \bf{R}^{m \times n}$，权重$λ$，$\varPhi_{0}$, $\varPhi_{-1} < -0, E_{0}, E_{-1}< -0, t_{0}, t_{-1}< -0$，$\mu_{0} \leftarrow 0.99\|\boldsymbol{D}\|_{2, 2}, \mu \leftarrow 10^{-5} \mu_{0}$。

当没有收敛时循环，否则退出循环。

$ \begin{array}{c} \boldsymbol{\varPhi}_{k} \leftarrow \boldsymbol{\varPhi}_{k}+\frac{\boldsymbol{t}_{k-1}}{t_{k}}\left(\boldsymbol{\varPhi}_{k}-\boldsymbol{\varPhi}_{k-1}\right)\\ \boldsymbol{E}_{k} \leftarrow \boldsymbol{E}_{k}+\frac{t_{k-1}}{t_{k}}\left(\boldsymbol{E}_{k}-\boldsymbol{E}_{k-1}\right)\\ \boldsymbol{Y}_{k}^{A} \leftarrow \widetilde{\boldsymbol{\varPhi}}_{k}-\frac{1}{2}\left(\widetilde{\boldsymbol{\varPhi}}_{k}+\widetilde{\boldsymbol{E}}_{k}-\boldsymbol{D}\right)\\ (\boldsymbol{U}, \boldsymbol{S}, \boldsymbol{V}) \leftarrow {svd}\left(\boldsymbol{Y}_{k}^{A}\right) \\ \boldsymbol{\varPhi}_{k+1} \leftarrow \boldsymbol{U}\left[\boldsymbol{S}-\frac{\mu}{2} \boldsymbol{I}\right]_+ \boldsymbol{V}^{*} \\ \boldsymbol{Y}_{k}^{E} \leftarrow \widetilde{\boldsymbol{E}}_{k}-\frac{1}{2}\left(\tilde{\boldsymbol{X}}_{k}+\widetilde{\boldsymbol{E}}_{k}-\boldsymbol{D}\right)\\ \boldsymbol{E}_{k+1} \leftarrow \operatorname{sign}\left[\boldsymbol{Y}_{k}^{E}\right] \times\left[\left|\boldsymbol{Y}_{k}^{E}\right|-\frac{\lambda \mu}{2} 11\right]_{+} \\ t_{k+1} \leftarrow \frac{1+\sqrt{1+4 t_{k}^{2}}}{2}, \mu \leftarrow \max (0.9 \mu, \bar{\mu}) \end{array} $

循环结束，输出$\boldsymbol{\varPhi, E}$。

3 实验

3.1 实验数据

训练样本集为从广州医科大学附属第一医院收集的多层螺线CT肺部扫描序列影像，共35例，其中2例的左肺壁上带有粘连病变。测试集来自EMPIRE10(Murphy等，2011)数据集，共22例，自带肺轮廓金标准，且有5例为肺壁上带有粘连病变。

图 1为训练样本集的部分影像，红色点为主标记点，绿色点为标记点。图 1(b)(c)中肺肿块分别长在肺区上叶中部和肺区下叶，可以看做是训练样本集中的噪声样本，因为它们肺轮廓线(绿色)走向与图 1(a)不同。

图 1 训练样本集影像的肺边缘及其标记点

Fig. 1 Training samples including noise samples are shown

((a) no juxta-pleural tumor; (b) with juxta-pleural tumor; (c) with juxta-pleural tumor)

3.2 实验步骤

1) 对训练集图像用灰度阈值法进行肺区分割，得到左肺边缘轮廓，将得到的肺轮廓的最上、最下、最左、最右及拐点作为主标记点。在每两个主标记点之间取19个点，共95个点作为标记点。如图 1所示。

2) 用RPCA训练形状模型，即对样本集标记点矩阵用式(4)、对协方差矩阵用式(8)进行低秩分解，并计算低秩分解后的协方差矩阵的特征值和特征向量，将特征值降序排列并选取与前$t$个最大特征值对应的特征向量，代入到形状模型$\boldsymbol{x}=\boldsymbol{\bar{x}}+\boldsymbol{P b}$中。基于含噪声样本的训练集分别用RPCA的ASM、PCA的ASM得到的形状模型，以及基于没有噪声样本的训练集，用PCA的ASM得到的形状模型分别如图 2—图 4所示，图中展示的是3种情况下前4个最大特征值对应的特征向量所生成的形状模型，蓝色轮廓为训练集样本的均值形状，红色轮廓为形状模型。

图 2 基于含噪声样本的训练集用RPCA的ASM得到的形状模型

Fig. 2 Shape models for training with noise samples by RPCA of ASM ((a) $t$ = 1; (b) $t$ = 2; (c) $t$ = 3; (d) $t$ = 4)

图 3 基于含噪声样本的训练集用PCA的ASM得到的形状模型

Fig. 3 Shape models for training with noise samples by PCA of ASM ((a) $t$ = 1; (b) $t$ = 2; (c) $t$ = 3; (d) $t$ = 4)

图 4 基于不含噪声样本的训练集用PCA的ASM得到的形状模型

Fig. 4 Shape models for training without noise samples by PCA of ASM ((a) $t$ = 1; (b) $t$ = 2; (c) $t$ = 3; (d) $t$ = 4)

3) 建立局部特征模型，本方法用的是梯度。对标记点$t$剖面邻域内的每个采样候选点$t$，都计算其梯度值$l_{k}$与标记点$t$的局部梯度模型中平均值$\bar{l}_{k}$之间的马氏距离

$f\left(l_{k}\right)=\left(l_{k}-\bar{l}_{k}\right)^{\mathrm{T}} \sum\limits_{k}^{-1}\left(l_{k}-\bar{l}_{k}\right)$

使$f(l_{k})$最小的点为标记点的下次迭代的位置，但在建立标记点的局部梯度模型时，需要剔除训练样本集中的噪声样本。具体的剔除算法为：求各个训练样本轮廓标记点的曲率，并用巴氏系数计算各自直方图的相似度，当相似度 < 0.8时，该直方图对应的样本为噪声样本。肺轮廓的曲率曲线如图 5所示。对于肺壁上无粘连肺肿块的肺轮廓标记点的曲率曲线如图 5(a)所示。对于肺壁上有粘连肺肿块的肺轮廓标记点的曲率曲线如图 5(b)所示，有无粘连肺肿块的曲率曲线图不同。

图 5 肺轮廓的曲率曲线

Fig. 5 Curvature curve of lung contour ((a) without juxta-pleural nodule; (b) with juxta-pleural nodule)

4 实验结果

4.1 所生成形状模型的比较

图 2—图 4中，当$t=1$时，表示仅用最大特征值对应特征向量来表达形状，可以看出图 2(a)与均值形状最近似，图 4(a)与均值形状较近似，图 3(a)与均值形状最不近似。当$t=2$时，表示用前两个特征值对应特征向量表达形状，可以看出图 2(b)和图 4(b)与均值形状都较近似，图 3(b)与均值形状最不近似。当$t=4$时，表示用前4个特征值对应特征向量表达形状，可以看出图 2(d)和图 4(d)已能很好地表达均值形状，但图 3(d)还是与均值形状不相符。

图 3中，图 3(a)—(c)中的红色形状已经不在正常的肺形状范围，而图 2和图 4中，从图 2(b)和图 4(b)开始，红色形状就已经在正常的肺形状范围，说明图 3(a)—(d)建立的形状模型受到了噪声样本的影响。而由于应用了RPCA，图 2(a)—(d)建立的形状模型没有受到噪声样本的影响，与无噪声样本的训练集建立的形状模型(图 4(a)—(d))基本相同。由此可以看出，RPCA能消除噪声样本对形状模型的影响。

4.2 分割结果评价

本文使用区域覆盖率(overlap)方法(孙申申等，2016)作为分割评价算法，具体为

$l_{p}=\frac{\boldsymbol{R}_{g} \cap \boldsymbol{R}_{n}}{\boldsymbol{R}_{g} \cup \boldsymbol{R}_{n}} \times 100 \%$

式中，$l_{p}$为区域覆盖率，$\boldsymbol{R}_{g}$和$\boldsymbol{R}_{n}$分别指金标准的像素集和分割算法的像素集。实验在训练样本集含有噪声样本的情况下，分割测试集的22个肺影像，具体结果为：基于RPCA的ASM的肺区分割的重叠度平均值为94.5%，而传统ASM(Sun等，2012)即基于PCA的ASM分割重叠度平均值为69.5%。与Shao等人(2014)的方法相比，由于也没有考虑训练集噪声样本的影响，故分割重叠度平均率为85%。将传统ASM(Sun等，2012)与基于RPCA的ASM分割效果相对比，如图 6所示。可以看出，基于RPCA的ASM分割方法能得到更好的分割结果。

图 6 ASM的肺区分割结果

Fig. 6 The lung segmentation result

((a) ASM based on PCA; (b) ASM based on RPCA)

图 7是用灰度阈值方法和基于RPCA的ASM分割方法分割含粘连肺壁肺肿块的影像结果。可以看出，基于RPCA的ASM分割方法能得到较好的分割结果。

图 7 含粘连肺壁型肺肿块的肺区分割结果

Fig. 7 Lung segmentation with juxta-peaural nodule

((a) traditional gray threshold segmentation method; (b) ASM based on RPCA method)

5 结论

本文主要研究了基于低秩理论的RPCA与ASM相结合的肺实质分割算法，比较了当训练样本集含有噪声样本时基于RPCA和基于PCA训练出的形状模型，发现基于RPCA的形状模型与不含噪声样本时基于PCA训练的形状模型相一致。基于RPCA的ASM的肺区分割的重叠度平均值为94.5%，而传统ASM即基于PCA的ASM分割重叠度平均值为69.5%。实验表明，基于PCA的传统ASM对非高斯型噪声敏感，基于RPCA的ASM分割肺区能得到比传统ASM更好的分割效果。

本文方法适用于正常的肺实质或粘连胸壁型肺实质分割。但对含大面积缺失的肺实质分割，本文方法并不适用，这是因为在边缘搜索过程中，异常标记点容易落到肺实质与病灶之间的过渡区域，而不是真正的肺边缘上，所以有必要识别异常标记点，并给异常标记点赋予不同的搜索函数，这是以后要解决的问题。

参考文献

Cheng X. 2013. Nonrigid Face Alignment for Unknown Subject in Video. Queensland: Queensland University of Technology

Cootes T F, Taylor C J, Cooper D H, Graham J. 1995. Active shape models-their training and application. Computer Vision and Image Understanding, 61(1): 38-59 [DOI:10.1006/cviu.1995.1004]

Hosseini-Asl E, Zurada J M, Gimel'farb G, El-Baz A. 2016. 3-D lung segmentation by incremental constrained nonnegative matrix factorization. IEEE Transactions on Biomedical Engineering, 63(5): 952-963 [DOI:10.1109/TBME.2015.2482387]

Murphy K, van Ginneken B, Reinhardt J M, Kabus S, Ding K, Deng X, Cao K L, Du K F, Christensen G E, Garcia V, Vercauteren T, Ayache N, Commowick O, Malandain G, Glocker B, Paragios N, Navab N, Gorbunova V, Sporring J, de Bruijne M, Han X, Heinrich M P, Schnabel J A, Jenkinson M, Lorenz C, Modat M, McClelland J R, Ourselin S, Muenzing A E A, Viergever M A. 2011. Evaluation of registration methods on thoracic CT:the EMPIRE10 challenge. IEEE Transactions on Medical Imaging, 30(11): 1901-1920 [DOI:10.1109/TMI.2011.2158349]

Oluyide O M, Tapamo J R, Viriri S. 2018. Automatic lung segmentation based on graph cut using a distance-constrained energy. IET Computer Vision, 12(5): 609-615 [DOI:10.1049/iet-cvi.2017.0226]

Peng Y G, Ganesh A, Wright J, Xu W L, Ma Y. 2012. RASL:robust alignment by sparse and low-rank decomposition for linearly correlated images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(11): 2233-2246 [DOI:10.1109/TPAMI.2011.282]

Qiu F, Hou F, Yuan Y, Wang W C. 2019. Blind image deblurring with reinforced use of edges. Journal of Image and Graphics, 24(6): 847-858 (邱枫, 侯飞, 袁野, 王文成. 2019. 加强边缘感知的盲去模糊算法. 中国图象图形学报, 24(6): 847-858D) [DOI:10.11834/jig.180543]

Shao Y Q, Gao Y Z, Guo Y R, Shi Y H, Yang X, Shen D G. 2014. Hierarchical lung field segmentation with joint shape and appearance sparse learning. IEEE Transactions on Medical Imaging, 33(9): 1761-1780 [DOI:10.1109/TMI.2014.2305691]

Sluimer I, Prokop M, van Ginneken B. 2005. Toward automated segmentation of the pathological lung in CT. IEEE Transactions on Medical Imaging, 24(8): 1025-1038 [DOI:10.1109/TMI.2005.851757]

Sofka M, Wetzl J, Birkbeck N, Zhang J D, Kohlberger T, Kaftan J, Declerck J and Zhou S K. 2011. Multi-stage learning for robust lung segmentation in challenging CT volumes//Proceedings of the 14th International Conference on Medical Image Computing and Computer-Assisted Intervention. Toronto: Springer: 667-674[DOI: 10.1007/978-3-642-23626-6_82]

Soliman A, Khalifa F, Elnakib A, El-Ghar M A, Dunlap N, Wang B, Gimel'farb G, Keynton R, El-Baz A. 2017. Accurate lungs segmentation on CT chest images by adaptive appearance-guided shape modeling. IEEE Transactions on Medical Imaging, 36(1): 263-276 [DOI:10.1109/TMI.2016.2606370]

Sun S H, Bauer C, Beichel R. 2012. Automated 3-D segmentation of lungs with lung cancer in CT data using a novel robust active shape model approach. IEEE Transactions on Medical Imaging, 31(2): 449-460 [DOI:10.1109/TMI.2011.2171357]

Sun S S, Fan L N, Kang Y, Ren H Z, Qi S L. 2016. Research on segmentation method of lung with juxta-pleural tumor based on the improved active shape model. Journal of Biomedical Engineering, 33(5): 879-884 (孙申申, 范立南, 康雁, 任会之, 齐守良. 2016. 基于改进主动形状模型的含胸壁粘连型肿块的肺区分割方法研究. 生物医学工程学杂志, 33(5): 879-884D) [DOI:10.7507/1001-5515.20160142]

van Ginneken B, Frangi A F, Staal J J, ter Haar Romeny B M, Viergever M A. 2002. Active shape model segmentation with optimal features. IEEE Transactions on Medical Imaging, 21(8): 924-933 [DOI:10.1109/TMI.2002.803121]

Wang G T, Zhang S T, Li F, Gu L X. 2013. A new segmentation framework based on sparse shape composition in liver surgery planning system. Medical Physics, 40(5): 051913 [DOI:10.1118/1.4802215]

Wright J, Peng Y G, Ma Y, Ganesh A and Rao S. 2009. Robust principal component analysis: exact recovery of corrupted low-rank matrices by convex optimization//Proceedings of the 22nd International Conference on Neural Information Processing Systems. Vancouver: Curran Associates Inc.: 2080-2088

Yang H, Cai S N, Wang J D and Quan L. 2014. Low-rank SIFT: an affine invariant feature for place recognition//Proceedings of 2014 IEEE International Conference on Image Processing. Paris: IEEE: 5731-5735[DOI: 10.1109/ICIP.2014.7026159]

Zhang S T, Zhan Y Q, Dewan M, Huang J Z, Metaxas D M, Zhou X S. 2012. Towards robust and effective shape modeling:sparse shape composition. Medical Image Analysis, 16(1): 265-277 [DOI:10.1016/j.media.2011.08.004]

Zhou X W, Yang C, Zhao H Y, Yu W C. 2014. Low-rank modeling and its applications in image analysis. ACM Computing Surveys, (CSUR), 47(2): 1-33 [DOI:10.1145/2674559]

Zhou X W and YU W C.2013. Low-rank modeling and its application in medical image analysis. Proc. SPIE 8750, Independent Component Analyses, Compressive Sampling, Wavelets, Neural Net, Biosystems, and Nanoengineering XI, 87500V, https://doi.org/10.1117/12.2017684