基于非负矩阵分解的多通道荧光显微图像光谱拆分方法与流程

1.本发明属于图像处理领域，更具体地，涉及基于非负矩阵分解的多通道荧光显微图像光谱拆分方法。


背景技术：

2.在多色成像实验中，每种荧光团通常会标记样品中的不同结构。在一张包含串扰的多通道数据叠加成的图像中，我们不知道某一像元中包含的信号究竟来自于哪一种荧光团，它可能就是某一种单独的荧光团发出的荧光，也有可能是几种荧光团所发出的荧光混合而得的结果，这取决于它们在靶细胞器或分子内的空间分布。光谱拆分就是将不同荧光团在光谱及空间上的重叠分离开来，确定每种荧光团对图像的每个像元的相对贡献，以反推出各种结构在样品上真实的位置及数量关系的一种数据处理方法。
3.目前被用于荧光显微成像中的光谱拆分算法多种多样，比较主流的方法分为两种思路：一种是基于分类的思想，它假设图像中所有像元均为纯像元，即图像中每一个像元的信号只来自于对应空间位置上的某一种荧光团所发出的荧光，表现在混合模型中即每个像元的丰度向量αi为一元向量(unary vector)，向量中只有一个元素为1，其余都为0；另一种则是基于参考光谱进行代数求解，假设图像中所有像元都是混合像元，图像中所有像元的灰度值都由样品中所包含的所有荧光团共同描述，但荧光团之间不存在相互作用，不同荧光团对于该像元的贡献可用丰度矩阵a来描述，与荧光团的激发以及发射特性有关。
4.目前所使用的光谱拆分都依赖于样品中所包含的荧光团光谱信息，想要实现多色的光谱拆分，常常需要在成像前经过比较繁复的先验实验针对荧光团获得较准确的光谱信息。考虑到样品制备周期以及荧光蛋白在生物体内的随机表达等问题，如果能跳过流程光谱测量这一步，简化光谱拆分的流程会是更好的选择。此外，在实际成像，尤其是对活体进行成像时，复杂的成像环境(如温度，酸碱度等)都会对荧光团的光谱产生影响，很难保证实际成像与光谱测量时环境条件绝对一致，所以拆分结果常常不理想。
5.在机器学习中，有一种称为无监督学习的方法，其输入数据时不需要附带标签即可进行分类，将这一思路转移到光谱拆分当中，直接从成像数据中计算光谱信息及荧光强度，不需要引入预先测量的光谱信息，可获取更优的光谱拆分结果。非负矩阵分解(non-negative matrix factorization，nmf)便是一种无监督学习的方法，无需使用任何先验信息(荧光团的激发谱或者发射谱)就能实现光谱拆分。不同于主成分分析(principal component analysis,pca)或是矢量分析(vector quantization,vq)，由于对矩阵的非负限制，数据在向量组合过程中只能使用加法而非减法，这样基于向量组合的分解具有可解释性和明确的物理意义，便称其分解出了局部特征。借助线性混合模型的矩阵表示形式p＝ea(e为端元矩阵，a为丰度矩阵，p即原图像)，nmf是依据最优化理论，通过迭代的方法，找到两个矩阵能够准确表述原图像，不需要提前测量光谱，在获取多通道图像之后将p进行矩阵分解，同时完成端元提取(单一荧光团图像)与丰度反演(荧光团光谱信息)。但是，由于丰度矩阵a与光谱信息有关，而不同荧光显微图像中的光谱信息会发生变化，因此丰度矩阵a的
稀疏度也会发生变化，但目前的nmf算法在迭代过程中并未考虑此因素。


技术实现要素：

6.针对现有技术的以上缺陷或改进需求，本发明提供了基于非负矩阵分解的多通道荧光显微图像光谱拆分方法，其目的在于解决现有的光谱拆分过程中图像存在串扰的技术问题。
7.为实现上述目的，按照本发明的一个方面，提供了基于非负矩阵分解的多通道荧光显微图像光谱拆分方法，包括以下步骤：
8.s1，获取多通道荧光显微图像，将每个通道的图像作为矩阵p中的列；
9.s2，将矩阵p分解为e与a两个矩阵，其中，e表示端元矩阵，a表示丰度矩阵；
10.s3，建立目标函数其中，f表示范数，α表示权重系数，j1(a)表示丰度矩阵a的稀疏度惩罚项，定义为j1(a)＝|sparseness(a)-spa|，sparseness(a)表示丰度矩阵a的稀疏度，spa表示丰度矩阵a的理想稀疏度值，0≤spa≤1；
11.s4，初始化端元矩阵e与丰度矩阵a；
12.s5，对目标函数中的端元矩阵e与丰度矩阵a进行迭代求解，直到目标函数小于预设值或迭代次数达到预设的最大迭代次数，则输出端元矩阵e和丰度矩阵a；
13.s6，将s5所得的端元矩阵e中的每列还原为荧光团图像，丰度矩阵a中的行表征为对应荧光团的光谱。
14.由于丰度矩阵a与光谱信息有关，当成像不同的荧光标记样本时，光谱信息会发生变化，因此丰度矩阵a的稀疏度也会发生变化，若直接令稀疏度趋近于零，那么当样品中存在共定位荧光标记时则容易得到像元过度分离的结果，导致图像失真严重。因此，通过上述技术方案，在目标函数上增加了丰度矩阵的稀疏度惩罚项，该稀疏度惩罚项中引入了丰度矩阵a的理想稀疏度值spa，并将其限定为[0，1]内的常数，通过减去spa，可以使丰度矩阵a的稀疏度值更接近理想值，保证丰度矩阵a是稀疏的，在迭代过程中实现精准约束，本技术中先使用荧光团现有光谱以及成像时所设置的荧光通道估算出每组实验中的丰度矩阵的理想稀疏度值，可在减少拆分图像串扰的同时，保证图像的连贯。
[0015]
优选地，s3中，0≤α≤m，m表示每个通道图像的像元数。
[0016]
优选地，s3中，l1表示一个矩阵的l1范数，定义为该矩阵中所有元素的绝对值之和，l2表示一个矩阵的l2范数，定义为该矩阵中所有元素的平方之和的平方根，n为丰度矩阵a的维数。
[0017]
优选地，s4中使用非负双奇异值分解对端元矩阵e与丰度矩阵a进行初始化。
[0018]
优选地，s4中初始化具体包括以下步骤：
[0019]
s401，对矩阵p进行部分奇异值分解，得到左奇异向量u、右奇异向量v以及对角矩阵s；
[0020]
s402，构建端元矩阵e的初始值e0，丰度矩阵a的初始值a0，j表示矩阵中的列：
[0021]
j＝1时，
[0022]
j≥2时，若||u
+
||||v
+
||＞||u-||||v-||，
[0023]
若||u
+
||||v
+
||＜||u-||||v-||，||，
[0024]
其中，t表示矩阵转置，uj,vj,sj分别表示左奇异向量u、右奇异向量v以及对角矩阵s的列，将左奇异向量u的列中所有≥0的元素保留，其余元素替换为0，记为u
+
，将左奇异向量u的列中所有《0的元素保留，其余元素替换为0，记为u-，将右奇异向量v的列中所有≥0的元素保留，其余元素替换为0，记为v
+
，将右奇异向量v的列中所有《0的元素保留，其余元素替换为0，记为v-。
[0025]
优选地，s5中迭代的具体步骤为：
[0026]
s501，将s402中的e0和a0作为迭代的初始值；
[0027]
s502，针对e，针对a，下降迭代步长μa：a:＝a-μae
t
(ea-p)；
[0028]
s503，将a的每一行投影至非负，保持其l2范数为1，将l1范数优化至spa；
[0029]
s504，判断目标函数是否小于预设值或迭代次数达到预设的最大迭代次数，达到则输出端元矩阵e和丰度矩阵a，否则重复s502、s503。
[0030]
优选地，还包括s7，采用光谱角距离和光谱信息散度对s6中获得的光谱进行相似度度量。
[0031]
优选地，还包括s8，根据s7中的相似度结果调节s3中的权重系数α和/或spa，采用s4、s5、s6获得新的荧光团图像和对应荧光团的光谱；
[0032]
s9，重复步骤s7-s8，直到s7中获得的相似度达到期望为止。
[0033]
优选地，s7中还包括，将s6中获得的荧光团图像和对应荧光团的光谱混合成反混图像，计算s1中的荧光显微图像和反混图像的均方根误差，以判断荧光显微图像光谱拆分后的信息损失量。
附图说明
[0034]
图1是基于非负矩阵分解的多通道荧光显微图像光谱拆分原理示意图；
[0035]
图2是一些实施例中基于非负矩阵分解的多通道荧光显微图像光谱拆分流程图；
[0036]
图3是又一些实施例中基于非负矩阵分解的多通道荧光显微图像光谱拆分流程图。
具体实施方式
[0037]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0038]
如图1和图2所示，本发明提出基于非负矩阵分解的多通道荧光显微图像光谱拆分方法，包括以下步骤：
[0039]
s1，获取多通道荧光显微图像，将每个通道的图像作为矩阵p中的列；
[0040]
s2，将矩阵p分解为e与a两个矩阵，其中，e表示端元矩阵，a表示丰度矩阵；
[0041]
s3，建立目标函数其中，f表示范数，α表示权重系数，j1(a)表示丰度矩阵a的稀疏度惩罚项，定义为j1(a)＝|sparseness(a)-spa|，sparseness(a)表示丰度矩阵a的稀疏度，spa表示丰度矩阵a的理想稀疏度值，0≤spa≤1；
[0042]
s4，初始化端元矩阵e与丰度矩阵a；
[0043]
s5，对目标函数中的端元矩阵e与丰度矩阵a进行迭代求解，直到目标函数小于预设值或迭代次数达到预设的最大迭代次数，则输出端元矩阵e和丰度矩阵a；
[0044]
s6，将s5所得的端元矩阵e中的每列还原为荧光团图像，丰度矩阵a中的行表征为对应荧光团的光谱。
[0045]
现有技术中的非负矩阵分解是以ea与p的均方根误差作为目标函数进行迭代优化，只考察了分解矩阵与原矩阵之间的相似性，因此很难在单次分解中将样品包含的所有荧光团同时拆分出来，且不同的图像之间会存在串扰，当存在串扰时，每种荧光团的信号依旧集中在少数通道之中，其余通道可视为0，这说明丰度矩阵a应该是稀疏的。而由于丰度矩阵a与光谱信息有关，当成像不同的荧光标记样本时，光谱信息会发生变化，因此丰度矩阵a的稀疏度也会发生变化。本技术通过在目标函数上增加了丰度矩阵的稀疏度惩罚项，并将该项定义为j1(a)＝|sparseness(a)-spa|，在常规的稀疏度的基础上引入了丰度矩阵a的理想稀疏度值spa，可根据查阅样品中所包含的荧光团光谱、荧光团激发效率以及多通道图像选用的滤光片进行估计，将其限定为[0，1]内的常数，通过减去spa，可以使丰度矩阵a的稀疏度值更接近理想值，既保证丰度矩阵a是稀疏的，在迭代过程中实现精准约束，最终减小了拆分后的荧光团图像之间的串扰，也避免了拆分图像中像元的过度分离，保留图像的连续性。
[0046]
具体地，将荧光显微图像做扁平化处理，然后将图像中的像元按照同样的顺序排成一维向量pi，i表示矩阵中的行，每个通道的图像组成矩阵p中的列。矩阵p分解即p＝ea，其中，e∈rm×k为端元矩阵，a∈rk×n为丰度矩阵，p∈rm×n即s1中的荧光显微图像，r表示矩阵中的所有元素均为实数，m为每个通道图像的像元数，n为图像采集通道数，k为样品中所包含的荧光团种类。根据s1和s2的矩阵分解形式建立目标函数，目标函数的第一项为ea与p的均方误差，当第二项加入本技术定义的稀疏度惩罚项后，可以根据不同的应用场景在0～1内选择spa的值，来保证丰度矩阵a为稀疏的，不受光谱信息变化的影响。
[0047]
进一步地，s3中，0≤α≤m，m为每个通道图像的像元数，在[0，m]内赋予稀疏度惩罚项一个权重系数，用于平衡分解相似性以及稀疏限制之间的权重，可以依据样品的荧光标记程度进行取值，当荧光样品本身为稀疏标记时，α的取值可接近m，而样品中本身荧光共定位的情况比较多时，α可靠近0取值。进一步，α的值可根据拆分评价后的结果进行调整。
[0048]
进一步地，s3中，l1表示一个矩阵的l1范数，定义为该矩阵中所有元素的绝对值之和，l2表示一个矩阵的l2范数，定义为该矩阵中所有元素的平方之和的平方根，n为丰度矩阵a的维数。因此，在本技术中丰度矩阵a的l1范数即
l2范数即当丰度矩阵a的稀疏度值为1时，表示a中只有一个非零元素，当值为0时，表示所有元素相等。
[0049]
进一步地，s4中使用非负双奇异值分解对端元矩阵e与丰度矩阵a进行初始化。具体包括以下步骤：
[0050]
s401，对矩阵p进行部分奇异值分解，得到左奇异向量u、右奇异向量v以及对角矩阵s；
[0051]
s402，构建端元矩阵e的初始值e0，丰度矩阵a的初始值a0，j表示矩阵中的列：
[0052]
j＝1时，
[0053]
j≥2时，若||u
+
||||v
+
||＞||u-||||v-||，||，
[0054]
若||u
+
||||v
+
||＜||u-||||v-||，||，
[0055]
其中，t表示矩阵转置，uj,vj,sj分别表示左奇异向量u、右奇异向量v以及对角矩阵s的列，将左奇异向量u的列中所有≥0的元素保留，其余元素替换为0，记为u
+
，将左奇异向量u的列中所有《0的元素保留，其余元素替换为0，记为u-，将右奇异向量v的列中所有≥0的元素保留，其余元素替换为0，记为v
+
，将右奇异向量v的列中所有《0的元素保留，其余元素替换为0，记为v
_
。
[0056]
目前在光谱拆分领域中，均通过随机产生非负矩阵的方式来初始化端元矩阵e与丰度矩阵a，这使得迭代常常只能优化值局部最优解，且每次分解得到的结果都不同，从而导致提取的单一荧光团图像的鲁棒性不强，且分解的荧光团光谱信息不准确。而本技术结合了非负双奇异值分解的思路，将其应用在初始化过程中，可以消除分解结果的随机性，也可以提升之后迭代收敛的速度。
[0057]
更具体地，基于非负双奇异值分解的思路，先对矩阵p进行首次部分奇异值分解，返回k个最大奇异值，得到左奇异向量um×k、右奇异向量vk×n以及对角矩阵sk×k，即pm×n≈um×ksk×
kvk
×
nt
，其中，对角矩阵中所有元素均为非负，该思路对应本技术的s401。
[0058]
因此，矩阵p可以写作向量形式p＝∑sjujvj·
或者p
(j)
＝s
j，j
×
uj×
vj·
。然后第二次奇异值分解的目标是在保留矩阵特征的情况下将矩阵中的元素均转化为非负，具体操作是将u
+
记为向量u中所有≥0的元素保留原值，其余元素用0替换，u-记为向量u中所有＜0的元素保留原值，其余元素用0替换，对向量v同理。随后，令x＝uv
t
＝(u
+-u
_
)(v
+-v-)
t
＝(u
+v+t
+u-v-t
)-(u
+
v-t
+u-v
+t
)＝x
+-x-，此时x中所有元素均为非负，那么p＝sx，p中所有元素均为非负，进一步，x可由x
+
近似：x≈x
+
。令μ
±
＝||u
±
||||v
±
||，那么x
+
可记为：即：当||u
+
||||v
+
||＞|||u-||||v-||，当||u
+
||||v
+
||＜|||u-||||v-||，该思路对应本技术的s402
[0059]
又由于p＝e0a0，在上述思路下对矩阵p进行分解，从而构建出：
[0060]
j＝1时：
[0061]
j≥2时：若||u
+
||||v
+
||＞||u-||||v-||，||，
[0062]
若||u
+
||||v
+
||＜||u-||||v-||，||，
[0063]
进一步地，在对端元矩阵e和丰度矩阵a初始化之后，开始进行迭代，s5中迭代的具体步骤为：
[0064]
s501，将s402中的e0和a0作为迭代的初始值；
[0065]
s502，针对e，针对a，下降迭代步长μa：a：＝a-μae
t
(ea-p)；
[0066]
s503，将a的每一行投影至非负，保持其l2范数为1，将l1范数优化至spa；
[0067]
s504，判断目标函数是否小于预设值或迭代次数达到预设的最大迭代次数，达到则输出端元矩阵e和丰度矩阵a，否则重复s502、s503。
[0068]
具体地，迭代过程中在s503利用了投影梯度下降的方法进行优化，可以在分解过程中精准约束稀疏度，结合本技术优化后的目标函数后，可以减少拆分后荧光图像之间的串扰。
[0069]
通过迭代收敛获得了最终输出的端元矩阵e与丰度矩阵a，然后在s6中还原成荧光团图像和对应的荧光团光谱。与其他的光谱拆分方法相比，本技术不需要提前假设多通道荧光图像数据种含有纯净端元，不需要提前测量样品中所包含荧光标记物的光谱数据，只需要设定端元数量便可直接进行光谱拆分；还可以同时获取多通道荧光图像的端元矩阵e与丰度矩阵a，不需要分步进行。
[0070]
由于本技术的技术方案无需提前获取光谱数据，因此，如图3所示，在另一些实施例中，可在s6之后进行拆分的精度评价，也就是进一步包括s7，采用光谱角距离和光谱信息散度对s6中获得的光谱进行相似度度量。
[0071]
尽管用于相似度度量的方法是现有的，但将其应用在本技术的s6之后，还可以根据该相似度结果来优化本技术s3中的目标函数，也就是进一步包括s8，根据s7中的相似度结果调节s3中的权重系数α和/或spa，采用s4、s5、s6获得新的荧光团图像和对应荧光团的光谱；s9，重复步骤s7-s8，直到s7中获得的相似度达到期望为止。通过对s6的拆分结果进行评价后，然后在α和spa的限定范围内调节目标函数中的α和/或spa的值，可以优化目标函数，直到获得期望达到的相似度为止。
[0072]
在s7中进行相似度度量的同时，还可以将s6中获得的荧光团图像和对应荧光团的光谱混合成反混图像，计算s1中的荧光显微图像和反混图像的均方根误差，以判断荧光显微图像光谱拆分后的信息损失量。该均方根误差越小，说明信息的损失量越少，由此可以进一步验证光谱拆分的精度。
[0073]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。