基于非负矩阵分解的显微图像彩色去卷积方法及分割方法

摘要

本发明公开了一种基于非负矩阵分解（NonnegativeMatrixFactorization，简称NMF）的显微图像彩色去卷积方法，属于图像信息处理技术领域。本发明方法针对染色标记后的组织显微图像，利用非负矩阵分解的方法将不同染色剂所对应的观测通道进行分离。本发明还公开了一种基于非负矩阵分解的显微图像分割方法，基于非负矩阵分解方法分离出的观测通道进行图像分割，不但可实现快速全自动化的图像处理，且获得的图像分割结果更准确，为后续的细胞检测、病理诊断分析提供了更准确的依据。相比现有技术，本发明能够有效提高显微图像分割的准确率，节省计算时间，且得到的显微图像的视觉效果更好。

说明书

技术领域

本发明涉及图像信息处理技术领域，尤其涉及一种基于非负矩阵分解(Nonnegative  Matrix Factorization，简称NMF)的显微图像彩色去卷积方法及分割方法。

背景技术

病理学是对细胞形态性质的微观研究。它在治疗方案的决策上起着主观重要的作 用。特别是在包括癌症在内的一些疾病的诊断中，病理图像的分析结果仍被认作为‘金 标准’。病理学的研究者已经认识到对病理图像定量分析的重要性。它可以用来支持临 床医生关于某种疾病的诊断决定，也可以帮助评价病人疾病的治疗效果。定量的分析不 仅对临床领域至关重要，同时对运用研究领域(如药品研发)也举足轻重。因此在病理 学中的计算机辅助分析能够加快临床医生的效率并提高他们诊断结果的准确性，最终造 福于广大的病人，减少他们的治疗代价甚至挽救他们的生命。病理切片组织样本或手术 样本是通过化学处理并分割到玻璃切片上的，其为疾病的诊治和预后治疗提供参考。数 字病理学则是对这些切片组织样本或手术样本在数字化的基础上进行的微观研究。为了 在显微镜下观测到多种组织成分，会使用几种不同的染色剂来对组织切片中的不同成分 进行染色。其中使用最广泛的为苏木精－伊红(H-E)染色，这种染色方法已经被病理学 家使用的100多年。苏木精把细胞核染成了深蓝色，而伊红把其他生物组织(如细胞质， 基质等)染成粉红色。细胞核由于其染色体和核仁的分布，容易显现出各种各样的图案。

与放射和其他形式的图像相比，处理病理组织显微图像的最大挑战在于其中巨大的 组织结构密度。举例来说，放大40倍的前列腺切片组织图像就有15000*15000个像素。 除了高密度，组织图像本身结构也很复杂，包含很多组织结构，而且这些结构没有明显 的边界，同一结构还存在很多的形态，这些都给处理带来巨大挑战。Wolberg针对乳腺 癌图像最早提出了基于主动轮廓模型的分割方法。在运用中它需要使用者手动输入细胞 核的位置，来方便程序进行分割。但是所遇到的问题是，位于基质但不在上皮层组织中 的细胞数量是十分巨大的。以前基于细胞或细胞核的检测是利用分割算法来实现的，这 些方法包括区域增长法、自适应阈值法、非监督的颜色聚类法，水平集法，形态学基础 的方法和基于监督的颜色纹理方法。但是这些方法对处理多变化形状和多尺度大小的人 体组织细胞时显得不是很鲁棒。这其中还存在组织切片获取，染色和数字化过程中人为 造成的不可避免的干扰现象。

对于组织切片图像来说，一般是用多种染色剂对其进行着色的，因为单一的着色会 使图像的部分信息丢失。在着色之后，如果能重新找到各染色剂所对应的观测通道，就 会使图像中的细胞信息更加清晰。因此，部分研究者会先利用彩色去卷积的方法将显微 图像中各染色剂所对应的观测通道分别提取出来，然后在此基础上进行图像分割，从而 获得更清晰准确的细胞或细胞核图像。例如，《基于着色分离的免疫组化图像核分割研 究》一文中提出一个基于着色分离的核分割方法,首先采用颜色去卷积算法对多着色的 IHC(immunohistochemical)图像进行着色分离,进而利用改进的SCFCM算法对单着色 灰度图像进行粗分割；然后利用分水岭算法分离粘连细胞；最后通过细胞核尺寸分析进行 后处理,完成对苏木素或多种染色的免疫组化图像的准确核分割。《Quantification of  histochemical staining by color deconvolution》这篇文章提出人工设计彩色去卷积，用于 分离苏木精，伊红和DAB三种染色。

然而，以上现有技术均利用人工设计彩色去卷积矩阵来实现观测通道的分离，需要 人为设计去卷积矩阵，而这个卷积矩阵的设计需要考虑染色成像时的光照强度以及染色 剂的成分信息，稍微的偏差就得不到理想的对应染色剂图像，因此难以实现图像处理的 实时性，且图像分割的准确率较低。

发明内容

本发明所要解决的技术问题在于克服现有技术通过人为介入的方式进行显微图像 彩色去卷积所带来的实时性差、准确率低的不足，提供一种基于非负矩阵分解的显微图 像彩色去卷积方法及一种基于非负矩阵分解的显微图像分割方法，能够有效提高显微图 像分割的准确率，节省计算时间，且得到的显微图像的视觉效果更好。

本发明具体采用以下技术方案解决上述技术问题：

基于非负矩阵分解的显微图像彩色去卷积方法，所述显微图像通过对经至多三种不 同染色剂染色标记的组织样本进行显微成像得到，该方法包括以下步骤：

步骤A、将大小为M×N的原始彩色显微图像的所有像素依次排列为行，以原始显微图 像中单个像素在RGB颜色空间中各分量所对应的光密度作为列，得到大小为 (M×N)×3的原始显微图像的光密度矩阵；

步骤B、以重构误差最小为目标，对所述光密度矩阵进行非负矩阵分解，将其分解为一 个(M×N)×K的矩阵A和一个K×3的矩阵S，K为所述染色剂的类别总数；

步骤C、根据下式得到第k种染色剂标记的光密度图像dAS_k，k＝1,2,…,K：

dAS_k＝dA(:,k)×dS(k,:)

式中，dA(:,k)表示矩阵A的第k列，dS(k,:)表示矩阵S的第k行；

步骤D、将第k种染色剂标记的光密度图像dAS_k在RGB颜色空间中重构，k＝1,2,…,K 得到第k种染色剂标记的光强度图像。

基于非负矩阵分解的显微图像分割方法，所述显微图像通过对经至多三种不同染色 剂染色标记的组织样本进行显微成像得到，首先利用如上所述彩色去卷积方法获取其中 至少一种染色剂标记的光强度图像，然后对所述光强度图像进行分割。

优选地，所述对光强度图像进行分割，包括利用阈值分割方法对光强度图像进行分 割的步骤。

进一步地，所述对光强度图像进行分割，还包括形态学处理的步骤，具体如下：对 阈值分割方法得到的分割图像依次进行以下形态学处理：膨胀、填充空洞、开运算-闭 运算、分裂。

相比现有技术，本发明方法具有以下有益效果：

(1)在同样的实验条件下，本发明方法的图像分割准确率比现有方法的准确率高；

(2)本发明寻找颜色去卷积矩阵是利用非负矩阵分解的方法自动实现，相比于现 有人工选择的方法，消除了人为因素的影响，更加可靠，快捷和方便，可实现完全自动 化；

(3)本发明使用形态学里的多种操作对得到的初步分割图像进行处理，最终获得 的分割结果更平滑，形状更精确。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的思路是针对染色标记后的组织显微图像，利用非负矩阵分解的方法将不同 染色剂所对应的观测通道进行分离，并基于分离出的观测通道进行图像分割，不但可实 现快速全自动化的图像处理，且获得的图像分割结果更准确，为后续的细胞检测、病理 诊断分析提供了更准确的依据。

下面以细胞图像分割为例来对本发明技术方案进行详细说明，其中组织切片样本采 用苏木精－伊红(H-E)染色，整个图像分割过程具体如下：

步骤1、构造新的二维图像矩阵：

将原使彩色显微图像的RGB三维矩阵(假设原图像大小为M×N个像素，则RGB 三维矩阵的大小为M×N×3)进行重新排列，将它构造为(M×N)×3的二维矩阵。行 为图像的像素总数，列为R、G、B三个通道的值，得到新的图像矩阵。

步骤2、对新的图像矩阵计算各像素R、G、B三个通道所对应的光密度，得到原始图 像的光密度矩阵；

在RGB模型中，R、G、B三个通道所对应的光强度I_R、I_G和I_B通过相机吸收后转 换为像素点，因为每个通道光强度的相关性是由一种在非线性方法下的着色浓度决定 的，图像的光强度值不可以直接用于分割和对每一种着色进行测量。R、G、B三个通 道的每个通道所对应的光密度(OD)可以表示如下：

${\mathrm{OD}}_C=-{\log }_{10}\left(\frac{I_C}{I_{0,C}}\right)=A_c*C_c---\left(\text{1}\right)$

其中，I_C是进入样本的光强度，I_0,C是观测到的光强度，A_c是对应某种吸收因子c 的染色强度，C_c是对应某种吸收因子c的R、G、B三个颜色的显示系数。可以看出， 对每个通道来说，OD和吸收的物质浓度是线性相关的，因此，在一个样本中，多种着 色的贡献可以用这个方法来分离。

步骤3、运用非负矩阵分解的方法分解图像的光密度矩阵，得到彩色去卷积矩阵：

现有技术均采用人工选择的方式来获取彩色去卷积矩阵，而本发明利用非负矩阵分 解的方法来自动获取，即以重构误差最小为目标，利用非负矩阵分解的方法将光密度矩 阵分解为一个(M×N)×K的矩阵和一个K×3的矩阵，K为所述染色剂的类别总数，本 实施例中K＝2，所得到的K×3的矩阵即为去卷积矩阵。

非负矩阵分解技术为现有技术，目前已在信号处理、模式识别、文本识别等领域得 到广泛应用。为了便于公众更好地理解，下面对其基本内容进行介绍。

在许多数据分析工作中，一个基本的问题是找到一个合适的数据代表，一个典型、 有用的数据代表可以使潜在的数据结构明确，且能够减少数据的维数。目前的研究中， 非负矩阵分解能够很好的解决这类问题。

对于一个给定矩阵X，矩阵A和矩阵S的最佳选择是要使得的X和AS之间的重 构误差最小，因此，NMF的求解问题实际上是个优化问题。

已有的非负矩阵分解算法大体可以分为三类：梯度下降算法、乘性迭代算法及交替 最小平方算法(ALS)。

给定A和S来观测X.考虑类似欧氏距离损失函数：

$J_E=D_E(X;\mathrm{AS})=\frac{1}{2}{\left|\right|X-\mathrm{AS}\left|\right|}_F^2=\frac{1}{2}\underset{\mathrm{pt}}{\Sigma }{(x_{\mathrm{pt}}-{\left[\mathrm{AS}\right]}_{\mathrm{pt}})}^2---\left(2\right)$

其中，X＝[x_pt]是样本中大小为p×t的观测量x_pt的数据阵， A＝[a_pn]＝[a₁,a₂,...,a_N]∈R^P×N是大小为p×n的混合矩阵。S＝[s_nt]是一个大小为n×t的 矩阵。对的一个简单的梯度下降步骤，我们希望通过下式对进行更新：

$S\leftarrow S-\eta \frac{{\partial J}_E}{\partial S}---\left(3\right)$

其中，η是一个很小的更新因子，上式也可以记作：

$s_{\mathrm{nt}}\leftarrow s_{\mathrm{nt}}-{\eta }_{\mathrm{nt}}\frac{{\partial J}_E}{{\partial s}_{\mathrm{nt}}}---\left(4\right)$

其中，对于(n,t)的每个结合，允许η_nt有不同的值。

为了计算局部的微分值，我们把目标函数记作：

$J_E=\frac{1}{2}{\left|\right|X-\mathrm{AS}\left|\right|}_F^2=\frac{1}{2}\mathrm{trace}\left({(X-\mathrm{AS})}^T(X-\mathrm{AS})\right)---\left(\text{5}\right)$

由于的一个无限小改变：

$S\leftarrow S+\partial S---\left(6\right)$

(2)中包含J_E的一个无限小改变：$J_E\leftarrow J_E+{\partial J}_E$

在(2)中，对这个无限小进行微分，得到：

$\left(\begin{array}{c}{\partial J}_E=-\mathrm{trace}({(X-\mathrm{AS})}^{T}A\partial S)\\ =-\mathrm{trace}(A^{T}X-A^T\mathrm{AS})\\ =-\underset{\mathrm{nt}}{\Sigma }{[A^{T}X-A^T\mathrm{AS}]}_{\mathrm{nt}}{\partial s}_{\mathrm{nt}}\end{array}\right)---\left(7\right)$

因此，

$\frac{{\partial J}_E}{{\partial s}_{\mathrm{nt}}}=-{[A^{T}X-A^T\mathrm{AS}]}_{\mathrm{nt}}=-({\left[A^{T}X\right]}_{\mathrm{nt}}-{\left[A^T\mathrm{AS}\right]}_{\mathrm{nt}})---\left(8\right)$

将(7)带入(3)，得到：

s_nt←s_nt+η_nt([A^TX]_nt-[A^TAS]_nt)  (9)

即s_nt＝[S]_nt的梯度更新步骤。由A和S的对称性，同理，得到a_pn＝[A]_pn的梯度更新步 骤：

a_pn←a_pn+η_pn([XS^T]_Pn-[ASS^T]_Pn)  (10)

一个简单的梯度更新算法在(9)和(10)之间交变直到收敛。为了维持a_pn和s_nt的所 有元素的非负性，我们使用下式表示：

s_nt←[s_nt+η_nt([A^TX]_nt-[A^TAS]_nt)]₊  (11)

其中，[S]₊＝max(0,s)是一个矫正方程，对a_pn也有同样的表示。

梯度更新只是一个简单的处理方法，其收敛很慢，且这种收敛性对每一步骤的尺寸 很敏感，为了克服这个问题，D.D.Lee和H.S.Seung提出了乘性迭代规则。

为了对s_nt构造一个乘性迭代规则，我们可以选择η_nt如：将(10)中的第一部分和 第三部分忽略掉，即s_nt＝η_nt[A^TAS]_nt或将它带入(10)，我们得到：

$s_{\mathrm{nt}}\leftarrow S_{\mathrm{nt}}\frac{{\left[A^{T}X\right]}_{\mathrm{nt}}}{{\left[A^T\mathrm{AS}\right]}_{\mathrm{nt}}}---\left(12\right)$

即的一个乘性更新规则。同理，可以得到的乘性迭代规则：

$a_{\mathrm{pn}}\leftarrow a_{\mathrm{pn}}\frac{{\left[{\mathrm{XS}}^T\right]}_{\mathrm{pn}}}{{\left[{\mathrm{ASS}}^T\right]}_{\mathrm{pn}}}---\left(13\right)$

与梯度下降算法相比，为了使(2)中，欧氏距离损失方程中的J_E更小，我们使用一 种类牛顿(Newton-like)方法找到交替的A和S，直接将J_E最小化。

对于一个固定矩阵A,首先，我们考虑对S的更新，将(8)用矩阵的形式表示，即：

$\frac{{\partial J}_E}{\partial S}=-(A^{T}X-A^T\mathrm{AS})---\left(14\right)$

它的最小值必须是0，则得到：

(A^TA)S＝A^TX  (15)

上式必须保持S使得J_E最小化，从(14)中,使用S＝-(A^TA)^-1A^TX，或用更加高效的线性 方程解法如MATLAB中的方程linsolve对S求解。同理，用A来最小化J_E时，我们可 以从方程(SS^T)A^T＝SX^T中解出A。

在ALS中，该算法不会自动给A和S强加非负约束，为了解决这个问题，最简单 的方法是在在正象限处理结果序列，从而得到最佳值:

S←[(A^TA)^-1A^TX]₊  (16)

A←[XS^T(SS^T)^-1]₊  (17)

其中，[M]₊将矩阵中所有的负值设置为0，由于在第一象限的处理，将负值移除，这意 味着负值将会没有理论意义，已经证明，这个过程在实际中有很好的应用。

步骤4、将分解后的光密度矩阵进行彩色去卷积算法重构，找到H通道对应的细胞图像；

本发明采用彩色去卷积算法找到细胞观测通道。各通道的图像重构按照以下公式：

dAS_H＝dA(:,1)×dS(1,:)  (18)

dAS_E＝dA(:,2)×dS(2,:)  (19)

其中，A表示非负矩阵分解的大小为(M×N)×2的矩阵；S表示2×3的系数矩阵，即彩 色去卷积矩阵；dAS_H表示H通道的光密度图像；dAS_E表示E通道的光密度图像，下 标H表示细胞图像通道，同理下标E表示细胞间质通道图像。

步骤5、强度重构；

在光密度图像基础上，分别重构H通道和E通道的光强度图像。根据Lambert-Beer 定律，通过对像素光密度的指数变换，可以得到光强度图像：

I_C＝exp(-d(AS))  (20)

其中，I_C是观测到的光强度，dAS表示光密度图像。

步骤6、细胞图像分割：

在步骤5中得到的细胞H通道光强度图像的基础上，进行细胞图像分割。此处可采 用现有的各种图像分割方法，例如期望最大化(EM)算法、分水岭算法、主动轮廓算 法等。本发明优选采用以下分割方法：

首先对H通道的光强度图像进行灰度化，并计算灰度化图像的平均灰度值；将此 灰度值作为阈值，对灰度图像进行阈值操作，高于阈值位置的像素值置为0，低于阈值 位置的像素值置为1，这样即得到初步的细胞分割结果。

在此时的分割图像中可以看到会有很多类似噪声的非细胞小块，可运用形态学中的 膨胀操作去除较小的噪声块。此时的结果里面只剩下细胞可能性较大的块，但还是可以 看到在这些大块中间还是存在许多盐粒的空洞，这时可采用形态学中的连通域内全填充 的方法来填补这些空洞，使得所有的块都是完整的、非缺失的。由于初步分割采用的是 阈值处理，所以图像中的可能细胞块边缘都比较尖锐，这时可采用形态学中的开运算- 闭运算来平滑这些区域。这时已经基本上得到细胞的分割结果，但考虑到图像中有很多 细胞块重叠的情况，可进一步采用在二值图像中寻找欧几里得局部最大值的方法来分裂 重叠细胞。经过分裂后的结果图即最终的细胞分割图，在最终分割图的基础上还可以标 记出细胞的位置。类似地，可利用E通道的光强度图像对伊红所标记的其他生物组织 图像进行分割和标记。

以上仅以两种染色剂着色的样本图像为例对本发明方法进行说明，实际上，对于任 意采用两种或三种染色剂着色的组织样本显微图像，采用本发明方法均可获得较好的分 割效果。

为了验证本发明方法的效果，进行了以下试验，并与现有的几种图像分割方法进行 了对比。

将本发明所提供的技术方案应用在用苏木精和苏木紫染色的乳腺癌组织彩色显微 图像上，图像的小为800×800。细胞分割过程具体如下：

1、重排像素矩阵：

原RGB图像矩阵大小为800×800×3，用reshape()函数重构为640000×3的新 的图像矩阵；

2、计算像素密度：

对得到的新的图像矩阵用(1)式进行计算各像素的光密度，得到原图像的光密 度矩阵；

3、对光密度矩阵进行非负矩阵分解：

对选取的图像的像素密度作为矩阵的行，使用(2)式中的代价函数来求解，直 到满足重构误差足够小；

4、颜色的分离：

采用彩色去卷积算法实现颜色的分离；分别找到对应的细胞图像-H通道和对应的 细胞间质密度图像-E通道。

5、各观测通道图像重构：

在密度图像基础上，用(20)式重构光强度图像。在新的光强度图像中，找到细胞 观测通道即H通道与细胞间质通道即E通道的光强度图像。

6、对H通道对应的光强度图像进行细胞图像分割：

(1)灰度化：

把H通道图像进行灰度化；

(2)阈值函数二值化

在灰度图的基础上利用阈值函数二值化图像，得到初步的分割结果：

(3)膨胀操作；

除去初步分割结果图中特别小的非细胞块；

(4)填补空洞；

由于采用的阈值分割，有些细胞块中会产生类似‘椒盐’的空洞噪声，这对最 后的分割效果会产生影响，所以利用连通域内无空缺的方法来填补这些空洞；

(5)开运算-闭运算；

消除块中狭窄的部分，去掉细的突出；

(6)把可能合并的细胞块分裂开；

在上面二值图像的基础上利用欧几里得距离求得局部最大值，在局部最大值的 基础上求得割图，即求得该分开的区域；根据分割的细胞图标记出细胞的位置。

为了验证本发明方法的效果，采用以下三种现有分割方法对上述实验原图进行细胞 分割、检测：对原图直接灰度化然后运用阈值分割；最大似然估计(EM)算法进行分 割；人工手动设计彩色去卷积矩阵进行分割，来与本发明方法进行比较，表1给出了4 种方法的实验结果：

表1 分割准确率对比

方法 DICE OL SN SP PPV 阈值分割 0.43 0.36 0.57 0.91 0.62 EM分割 0.62 0.53 0.59 0.94 0.73 手动设计参数彩色去卷积 0.65 0.52 0.62 0.96 0.79 本发明方法 0.77 0.63 0.69 0.99 0.89

从表1中可以看出，采用同样的后续分割方法时，本发明方法得到的分割准确率高 于其他几种模型的分割结果。

此外，根据实验所得的分割图像可以看出，本发明方法得到的细胞通道图像比手动 设计参数的彩色去卷积矩阵得到的图像更能突出原图中细胞的信息。

综上可知，相对于现有的图像分割方法，本发明方法无论是从准确率还是从可实行 效果上来看，都具有明显的优势。