CN103236052B - 一种基于l1最小化模型的全自动细胞定位方法 - Google Patents

Abstract

本发明公开了一种基于L1最小化模型的全自动细胞定位方法，步骤为：①对原始图像进行二值化处理，得到二值化图像B；②对B进行腐蚀，得到包括神经元的连通域，并对每个连通域按③和④进行处理，得到所有细胞的位置；③对每个连通域，将其嵌入到一个长方体，称为子区域，从子区域中获取候选种子点，再根据候选种子点间的距离剔除冗余的候选种子点，将剩余的候选种子点作为真实种子点，真实种子点的个数记为k；④构建L1最小化模型，得到连通域中真实的细胞位置和半径。本发明采用提取子区域的方式进行细胞定位，每个子区域可同时独立进行运算，可以将所有连通域并行运算，可极大提高整体运算速度，因而能有效处理大范围数据。

Description

一种基于L1最小化模型的全自动细胞定位方法

技术领域

本发明属于生物医学光学图像处理领域，特别是涉及生物医学光学图像中细胞定位方法。具体为一种基于L1最小化模型的全自动细胞定位方法，其中，L1为1范数。本发明方法尤其适用于神经元胞体的定位。

背景技术

神经回路是大脑功能的物理基础。绘制精细神经回路图能极大地提高人类关于大脑功能的认识水平。如果我们预先完成定位神经细胞的工作，那么我们将能够更快更精确地追踪神经突起，从而高效完成神经回路图的绘制工作。事实上这一方法已经在神经科学界得到广泛的应用，例如用这种方法我们成功研究了癌症干细胞与神经微环境的依赖关系，成功找到了神经细胞与神经血管两者之间的分布规律。最近突飞猛进的荧光标记和显微成像技术，使得人类能够以微米甚至超微米量级的水平测量如老鼠之类的啮齿类动物的完整大脑回路结构。不过即使是绘制仅仅像老鼠这样的小动物的大脑神经回路图，人类就要处理浩瀚如海的数据，这是人类大规模分析大脑数据所面临的巨大障碍。

自动定位和分割三维神经图像细胞的方法取得了巨大进展，典型的方法有分水岭算法，梯度流算法，这些方法只适用于简单形态的细胞的定位和分割。最近FARSIGHT的方法【参见文献：1.Al-Kofahi,Y.,Lassoued,W.,Lee,W.&Roysam,B.Improved automatic detection and segmentationof cell nuclei in histopathology images.IEEE Trans.Biomed.Eng.57,841-852(2010).】引人注目，它是一种多尺度滤波的方法，它利用不同尺度下的模板与图像匹配，达到定位胞体的目的。然而，该方法对于具有粗突起的神经元细胞不能很好地处理。神经细胞的定位工作仍然要依赖人工辅助，计算机自动定位神经细胞的方法还远远不够成熟。

发明内容

本发明的目的是提供一种基于L1最小化模型的全自动细胞定位方法，该方法处理速度快，并且可以在粗突起的干扰下实现细胞尤其是神经元胞体定位。

本发明提供的一种基于L1最小化模型的全自动细胞定位方法，包括下述步骤：

第1步、读取原始图像，并进行二值化处理，得到二值化图像B；

第2步、对二值化图像B进行腐蚀，得到包括神经元细胞的连通域，组成连通域集；对该连通域集中的每一个连通域均按照下述第3步与第4步的过程处理，得到所有细胞的位置和半径；

第3步、将提取的连通域嵌入到一个长方体区域，称之为子区域，从该子区域中获取候选种子点，再根据候选种子点间的距离剔除冗余的候选种子点，将剩余的候选种子点作为真实种子点，真实种子点的个数记为k；

第4步、构建如式I所示的L1最小化模型，得到连通域中真实的细胞位置和半径；其具体过程为：

第4.1步、求解式I中的L1最小化模型以估计参数，得到各真实种子点所对应的位置和半径：

$\min {\left|\right|\underset{o\∈V}{\mathrm{\Σ}}{B}_L\left(o\right)-\underset{o\∈V}{\mathrm{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}f(o,o_j,r_j)\left|\right|}^{1/3}+\mathrm{\λ}{\mathrm{\Σ}}_{j=1}^k{r}_j$    式I

s.t.r_j≥0,j＝1,2,…,k

B_L是指所述子区域的二值化图像，V为B_L的坐标集，o为B_L上的像素点，o_j为第j个真实种子点的坐标值，r_j是为第j个真实种子点的半径；f(o,o_j,r_j)表示以o_j为球心r_j为半径的球函数，λ为式I中前、后两部分的权衡，即权值，0＜λ＜1；

第4.2步、对于所有真实种子点所对应的半径，如果其大于等于原始图像中最小细胞的半径，则认为此半径对应的位置为真实的细胞位置，保存此细胞位置和半径。

作为上述技术方案一种改进，第1步具体包括下述过程：

第1.1步、读入三维原始图像，记为I；

第1.2步、将I与预设的阈值T1做取小操作，然后将之与一个均值模板做多次卷积操作，得到一副新图像，记为C；T1的取值为0.8—0.9倍图像I的背景的平均像素值；

第1.3步、根据如下核函数，利用I和C，得到二值化图像B：

其中，o表示三维图像上的像素点，I(o)和C(o)分别表示I和C在点o处的像素值。

作为上述技术方案另一种改进，第2步具体包括下述过程：

第2.1步、令R=B，对于R上的每一个像素点，如果该点与其周围26个点的二值化像素值之和小于阈值T2，则将此点的值置为0，否则其值保持不变，将进行上述操作后的图像记为E；T2的取值范围为8～14个像素；

第2.2步、对于图像E，提取其中像素个数在预设的范围内之间的连通域，并存储在连通域集待用；该预设的范围根据最小细胞与最大细胞所占像素个数确定；

第2.3步、将图像E中已经提取出来的连通域的像素值置为0，得到新的图像，用之来替换R；

第2.4步、对于新得到的R，重复第2.1步～第2.3步直至连通域的体素个数小于预设的范围的下限。

作为上述技术方案再一种改进，第3步具体过程为：

第3.1步、从连通域集中提取一个连通域，将提取出来的连通域嵌入到一个长方体区域，将之称为子区域，并将该子区域的二值化图像和原始图像分别记为B_L和I_L；

第3.2步、将B_L和I_L同时与全1模板进行卷积，得到的图像分别记为B⁺和I⁺；

第3.3步、根据图像B⁺和I⁺，从子区域中获取候选种子点；

第3.4步、从候选种子点中根据候选种子点间的距离剔除冗余的候选种子点，将剩余的候选种子点作为真实种子点。

作为上述技术方案进一步改进，第2.2步中，所述预设的范围为L_min～L_max，L_min=100±20像素，L_max=20000±200像素。

作为上述技术方案进一步改进，第3.3步中，对于子区域里面的每一个点o，如果该点的像素值B⁺(o)大于预设的阈值T3，而且I⁺(o)为局部极大值，将此点作为候选种子点，T3的取值范围为60～100个像素。

作为上述技术方案进一步改进，第3.4步的具体过程为，从候选种子点中选出I⁺(o)最大的那个点，将该点记为o₁，剔除与点o₁的欧拉距离小于预设的阈值T4的候选种子点，重复此操作直至候选种子点集为空集或重复至少10次，将剩余的候选种子点作为真实种子点，T4的取值为最小细胞的半径的0.8～0.9倍。

作为上述技术方案进一步改进，其特征在于，第4.1步中的式I转化为式II：

$\min {\left|\right|\underset{o\∈V}{\mathrm{\Σ}}{B}_L\left(o\right)-\underset{o\∈V}{\mathrm{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}f(o,o_j,r_j)\left|\right|}^{1/3}+\mathrm{\λ}{\mathrm{\Σ}}_{j=1}^k{w_{j}r}_j$    式II

s.t.r_j≥0，j＝1，2，…，k

式中，w_j为权重。

作为上述技术方案更进一步改进，对于式II，用梯度投影算法来计算半径r_j，利用均值算法来计算位置o_j。

上述技术方案尤其适用于神经元胞体的定位。

基于L1最小化模型的全自动细胞定位方法与现有技术相比，本发明具有以下技术特点：

1、能在粗突起的干扰下实现细胞尤其是神经元胞体的定位。

本发明方法不仅能够对各类细胞进行定位，而且特别适用于有粗突起的神经元细胞的定位。

对于每一个提取的子区域，本发明首先获取多个种子点，作为神经元胞体位置的预估，由于神经元在三维空间中的稀疏性，所以很多预估位置是错误的，这些错误的位置所对应的球半径应该趋于0，这一点也体现出了胞体半径序列的稀疏性，本发明采用L1最小化模型表征半径序列的稀疏性。本发明利用球函数对子区域进行拟合，同时考虑半径序列稀疏性，构成一个基于L1最小化模型，通过对该模型的求解，就可以使得错误的预估位置半径趋于0，这样就可以在粗突起的干扰下较准确地实现神经元胞体的定位。

2、本发明方法提出了一种L1最小化模型，该模型可以联合利用梯度投影算法与均值算法进行求解，有效地提升了计算速度，以实现快速定位。

对于需要求解的L1最小化模型，每个种子点都会对应4个参数（3个位置参数和1个半径参数），所以参数很多，如果将所有的参数作为一个整体，用梯度投影算法来进行求解，将会十分耗时；本发明将参数分为位置和半径两个部分，用梯度投影算法来计算半径，用均值算法来计算位置，由于均值算法的计算速度比梯度投影算法的计算速度要快得多，这样就能有效的提升算法的速度。

3、本发明方法能有效处理大范围数据。

本发明采用提取子区域的方式来进行细胞定位，每个子区域可以同时独立进行处理，因此可以将所有连通域并行处理，可极大提高整体处理速度，因而能有效处理大范围数据。

本发明方法为计算机全自动神经细胞定位和全脑神经突起追踪与重建提供重要前提条件。

附图说明

图1为本发明实例提供的全自动胞体定位方法的流程图；

图2是本发明实例提供的三维原始图像；

图3是本发明实例中三维原始图像合并像素点后得到的图像；

图4是本发明实例中合并像素点后图像的背景图像；

图5是本发明实例中合并像素点后图像的二值化图像；

图6是本发明实例中二值化图像的连通域集图像；

图7是本发明实例中二值化图像的残余图像；

图8是本发明实例提取的某一个连通域所对应的子区域的图像；

图9是本发明实例中子区域所获得的候选种子点；

图10是本发明实例中子区域所获得的真实种子点；

图11是本发明实例中子区域所确定的真实胞体位置；

图12是本发明实例中原始图像的所有的真实胞体的定位位置。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明。在此需要说明的是，对于这些实施方式的说明用于帮助理解本发明，但并不构成对本发明的限定。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示，本发明具体实施方式提供的一种全自动细胞定位方法的具体过程为：

第1步、读取神经元胞体的原始图像，并进行二值化处理，得到二值化图像B；

本发明对获取原始图像的二值化图像方式并没有特别要求，下面只是示例一种方式，其具体过程为：

第1.1步读入三维原始图像，记为P；

第1.2步将P与一个阈值T1(T1为0.8～0.9倍图像I的背景的平均像素值)做取小操作，然后将之与一个均值模板做多次（如10-30次）卷积操作，得到一副新图像，记为C；

第1.3步根据如下核函数（式①），利用P和C，可以得到二值化图像，记为B：

其中，o表示三维图像上的像素点，P(o)和C(o)分别表示P和C在点o处的像素值。

第2步对二值化图像B进行腐蚀，得到包括神经元的连通域，组成连通域集；

由于神经元突起分布的密集性和形态的复杂性，通常会形成很大的连通域，这将使计算变得复杂且损耗大量计算资源及时间，为此，我们需要将二值化图像B进行腐蚀，将大的连通域腐蚀成小的连通域。下面示例其中一种方式：

第2.1步在进行腐蚀操作前，先给出初始腐蚀条件，令图像R=B；

对于R上的每一个像素点，如果该点与其周围26个点的二值化像素值之和小于阈值T2(T2的取值范围为8～14)，则将此点的值置为0，否则其值保持不变，将进行上述操作后的图像记为E；

第2.2步对于E，提取其中像素个数在L_min～L_max之间的连通域，并存储在连通域集待用；

L_min～L_max取值范围根据最小神经元与最大神经元所占像素个数来估计，L_min的取值范围为100±20，L_max的取值范围为20000±200；

第2.3步将E中已经提取出来的连通域的像素值置为0，得到新的图像，用之来替换R，对于新的R；

第2.4步重复第2.1步～第2.3步直至连通域的体素个数小于L_min。

对连通域集中的每一个连通域均按照下述第3步与第4步的过程处理，得到所有胞体的位置；

第3步从连通域集中提取一个连通域，将提取出来的连通域嵌入到一个长方体区域，将之称为子区域，从该子区域中获取候选种子点，再根据候选种子点间的距离剔除冗余的候选种子点，将剩余的候选种子点作为真实种子点，真实种子点的个数记为k。

提取出来的连通域中，可能包含0个、1个或则多个真实的胞体。在寻找胞体体位置之前，先获取种子点，作为胞体位置的预估。下面示例说明其具体实现过程：

第3.1步从连通域集中提取一个连通域，为了方便后续处理，我们将提取出来的连通域嵌入到一个长方体区域，将之称为子区域，并将该子区域的二值化图像和原始图像分别记为B_L和I_L。

第3.2步将B_L和I_L同时与一个全1模板进行卷积，得到的图像分别记为B⁺和I⁺。

第3.3步对于子区域里面的每一个点o，如果该点的像素值B⁺(o)大于阈值T3而且I⁺(o)为局部极大值，我们就将此点作为候选种子点；

T3由全1模板像素点个数决定，此处取值范围为60～100；

第3.4步为了计算简单，我们需要删除一些冗余的种子点，从候选种子点中选出I⁺(o)最大的那个点，将该点记为o₁，剔除与点o₁的欧拉距离小于阈值T4的候选种子点，重复此操作直至候选种子点集为空集或重复至少10次，将剩余的候选种子点作为真实种子点，真实种子点的个数记为k；

T4的取值应接近最小胞体的半径，如三维原始图像中最小胞体的半径的0.8～0.9倍；

第4步、构建L1最小化模型，得到一个连通域中胞体的位置；

此时已经获得了此子区域的真实种子点（胞体位置的预估），本方法将用第4.1步～第4.3步来剔除错误的预估位置，获得真实的胞体位置；

第4.1步将B_L用以每个真实种子点为球心，给定的初值为半径的球函数来做拟合，可得如下表达式②

$B_L\left(o\right)={\mathrm{\Σ}}_{j=1}^{k}f(o,o_j,r_j)+\mathrm{res}\left(o\right)$    式②其中，f(o,o_j,r_j)表示第j个球函数，其表达式为式③：

$f(o,o_j,r_j)=\left\{\begin{array}{cc}1& \mathrm{if}|o-o_j|\≤r_j\\ \exp (-\frac{{\left(\right|o-o_j|-r_j)}^2}{0.1})& \mathrm{otherwise}\end{array}\right.$    式③

res(o)为用所有球函数拟合B_L后没有被拟合到的部分，o为B_L上的像素点，k为真实种子点的个数，o_j为第j个种子点的坐标值，r_j为第j个种子点的半径；

由于B_L包含了粗突起的部分，前面获得的k个胞体的预估（种子点)位置(o₁，o₂，…，o_k)中，有一部分是错误的胞体位置，所以神经元胞体定位这一问题就转化为了从(o₁，o₂，…，o_k)中剔除错误的位置，寻找真实的胞体位置；

以往通常会采用非线性最小二乘法来求解，也就是说通过求解如下最小剩余函数式④这一优化问题来估计参数

$\min \left|\right|\underset{o\∈V}{\mathrm{\Σ}}{B}_L\left(o\right)-\underset{o\∈V}{\mathrm{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}f(o,o_j,r_j)\left|\right|$    式④

s.t.r_j≥0,j＝1,2,…,k

其中，||||为2-norm，V为B_L的坐标集；

然而在有粗突起的干扰下，通过求解式④得到的胞体位置并不理想，为了解决这一问题，本发明根据神经元的稀疏性，我们提出用L1最小化模型来解决这一问题，模型表达式如下式⑤：

$\min {\left|\right|\underset{o\∈V}{\mathrm{\Σ}}{B}_L\left(o\right)-\underset{o\∈V}{\mathrm{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}f(o,o_j,r_j)\left|\right|}^{1/3}+\mathrm{\λ}{\mathrm{\Σ}}_{j=1}^k{r}_j$    式⑤

s.t.r_j≥0，j＝1，2，…，k

λ为式⑤前、后两部分的权衡(取值范围为0.01～0.5，如取为0.025)。

第4.2步求解优化问题式⑤等同于求解如下优化问题式⑥：

$\min {\left|\right|\underset{o\∈V}{\mathrm{\Σ}}{B}_L\left(o\right)-\underset{o\∈V}{\mathrm{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}f(o,o_j,r_j)\left|\right|}^{1/3}+\mathrm{\λ}{\mathrm{\Σ}}_{j=1}^k{r}_j$   式⑥

s.t.r_j≥0，j＝1，2，…，k

求解优化问题式⑥可以得到k个位置和半径，其求解过程如下：S1设置初始权重w_j＝1，j＝1，2，…，k，设定r_j的初始值为2～10个像素；S2令t＝0，t表示迭代次数，设初值

)其中

为1,2,…,k个种子点的坐标值，

为第1，2，…，k个种子点设定的初始半径；

S3根据用梯度投影算法获得

所述梯度投影算法的步骤如下：

为方便说明，记原来式⑦：

$\min {\left|\right|\underset{o\∈V}{\mathrm{\Σ}}{B}_L\left(o\right)-\underset{o\∈V}{\mathrm{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}f(o,o_j,r_j)\left|\right|}^{1/3}+\mathrm{\λ}{\mathrm{\Σ}}_{j=1}^k{w_{j}r}_j$    式⑦

S3.1取初始可行点

m表示迭代次数；

S3.2利用式⑧计算梯度

$g_i^m=\left\{\begin{array}{cc}\&dtri;h\left(r_i^m\right)& \mathrm{if}{r}_i^m>0\mathrm{or}\&dtri;h\left(r_i^m\right)<0\\ 0& \mathrm{otherwise}\end{array}\right.$    式⑧S3.3计算初始步长 ${\mathrm{\&alpha;}}_0=\underset{\mathrm{\&alpha;}}{\arg \min \mathrm{}}h(r^m-{\mathrm{\&alpha;g}}^m)$

在序列α₀,βα₀,β²α₀,…,β＝1/2选取第一个满足下式的值，记为α^m

$h(r^m-{\mathrm{\&alpha;}}^m\&dtri;h\left(r^m\right))\&le;h\left(r^m\right)-\mathrm{\&mu;}\&dtri;f{\left(r^m\right)}^T({\mathrm{\&alpha;}}^m\&dtri;h\left(r^m\right))$

S3.4令 $r^{m+1}=r^m-{\mathrm{\&alpha;}}^m\&dtri;h\left(r^m\right)$

S3.5如果收敛，则停止；否则令m＝m+1，转入S3.1。

S4根据

用均值法获得

具体计算如式⑨：

对于每一组参数

计算

   式⑨

其中，以

为球心，以

为半径的球的内部区域的像素点的个数记为n_j，其中像素点的坐标记为

是像素点

坐标的平均值。

S5如果结果收敛（收敛条件为所有参数前后两次迭代之间的变换小于0.001），则转入S6，否则令t+1→t，重复S3～S5

S6利用式⑩更新权重

$w_j=\min \{\frac{\max (r_1,r_2,...,r_k)}{r_j+0.01},40\},j=\mathrm{1,2},...,k$    式⑩

S7重复S2～S6直至所有的参数收敛。

第4.3步通过上述对优化问题式⑥的求解，可以得到k个位置和半径，对于所有真实种子点所对应的半径，如果其大于等于原始图像中最小胞体的半径的，则认为此半径对应的位置为真实的胞体位置，保存此胞体位置和半径。

实例：

以超分辨荧光成像显微镜或功能性双光子共聚焦成像显微镜获取的老鼠大脑切片图像为对象，对原始图像进行预处理，以方便后续的各项操作。

第1步：

读入三维原始图像（如图2所示），将三维原始图像的每一帧二维图像的4个(每个田字格)像素合并成一个像素，每个像素的信号值采用直接相加的方法，将得到的新的图像记为I（如图3所示）；

将I与T1(T1取400为宜)做取小操作，然后将之与一个9x9x1的均值模板做20次卷积操作，得到的新图像称为背景图像，记为C（如图4所示）；

根据式Ⅰ，利用I和C可以得到二值化图像，记为B（如图5所示）；

第2步：

按照上述具体实施方式中第2步的流程，将其中各参数按如下取值，T2＝27×(1/3+0.027×i)，L_min=100，L_max=20000，所得到的连通域集和残余图像分别如图5和图6所示；

第3步：

第3.1步从连通域集中提取一个连通域，为了方便后续处理，我们将提取出来的连通域嵌入到一个长方体区域，将之称为子区域(如图8所示)，并将该子区域的二值化图像和原始图像分别记为B_L和I_L。

第3.2步将B_L和I_L同时与一个7x7x5的全1模板进行卷积，得到的图像分别记为B⁺和I⁺。

第3.3步取T3=70，得到的候选种子点如图9所示；

第3.4步取T4=4.8μm，重复次数为10，得到的真实种子点如图10所示；

第4步：

第4.1步：第j个种子点设定的初始半径为5个像素，λ取为0.025，通过式Ⅵ的求解，得到k个位置和半径；

第4.2步T5=2.5，得到的真实的胞体位置和半径，如图11所示；

第5步：

所得到的全部的真实的胞体位置，如图12所示。

上述内容是以神经元胞体为例进行说明的，本发明可以用于处理带突起与不带突起的各类神经元胞体。由于普通细胞的形状与不带突起的神经元胞体是类似的，故本发明所提供的技术方案同样适用于其它各种细胞的定位。

以上所述为本发明的较佳实施例而已，但本发明不应该局限于该实施例和附图所公开的内容。所以凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。

Claims (10)

1.一种基于L1最小化模型的全自动细胞定位方法，包括下述步骤：

第1步、读取原始图像，并进行二值化处理，得到二值化图像B；

第2步、对二值化图像B进行腐蚀，得到包括细胞的连通域，组成连通域集；对该连通域集中的每一个连通域均按照下述第3步与第4步的过程处理，得到所有细胞的位置和半径；

第3步、将提取的连通域嵌入到一个长方体区域，称之为子区域，从该子区域中获取候选种子点，再根据候选种子点间的距离剔除冗余的候选种子点，将剩余的候选种子点作为真实种子点，真实种子点的个数记为k；

第4步、构建如式I所示的L1最小化模型，得到连通域中真实的细胞位置和半径；其具体过程为：

第4.1步、求解式I中的L1最小化模型以估计参数，得到各真实种子点所对应的位置和半径：

$\min {\left|\right|\underset{o\&Element;V}{\mathrm{\&Sigma;}}{B}_L\left(o\right)-\underset{o\&Element;V}{\mathrm{\&Sigma;}}\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}f(o,o_j,r_j)\left|\right|}^{1/3}+\mathrm{\&lambda;}{\mathrm{\&Sigma;}}_{j=1}^k{r}_j$    式I

s.t.r_j≥0,j＝1,2,…,k

B_L是指所述子区域的二值化图像，V为B_L的坐标集，o为B_L上的像素点，o_j为第j个真实种子点的坐标值，r_j是为第j个真实种子点的半径；f(o,o_j,r_j)表示以o_j为球心r_j为半径的球函数，λ为权值，0＜λ＜1；

第4.2步、对于所有真实种子点所对应的半径，如果其大于等于原始图像中最小细胞的半径，则认为此半径对应的位置为真实的细胞位置，保存此细胞位置和半径。

2.根据权利要求1所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，第1步具体包括下述过程：

第1.1步、读入三维原始图像，记为P；

第1.2步、将P与预设的阈值T1做取小操作，然后将之与一个均值模板做多次卷积操作，得到一副新图像，记为C；T1的取值为0.8—0.9倍图像I的背景的平均像素值；

第1.3步、根据如下核函数，利用P和C，得到二值化图像B：

其中，o表示三维图像上的像素点，P(o)和C(o)分别表示P和C在点o处的像素值。

3.根据权利要求1所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，第2步具体包括下述过程：

第2.1步、令R=B，对于R上的每一个像素点，如果该点与其周围26个点的二值化像素值之和小于阈值T2，则将此点的值置为0，否则其值保持不变，将进行上述操作后的图像记为E；T2的取值范围为8～14个像素；

第2.2步、对于图像E，提取其中像素个数在预设的范围内之间的连通域，并存储在连通域集待用；该预设的范围根据最小细胞与最大细胞所占像素个数确定；

第2.3步、将图像E中已经提取出来的连通域的像素值置为0，得到新的图像，用之来替换R；

第2.4步、对于新得到的R，重复第2.1步～第2.3步直至连通域的体素个数小于预设的范围的下限。

4.根据权利要求1所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，第3步具体过程为：

第3.1步、从连通域集中提取一个连通域，将提取出来的连通域嵌入到一个长方体区域，将之称为子区域，并将该子区域的二值化图像和原始图像分别记为B_L和I_L；

第3.2步、将B_L和I_L同时与全1模板进行卷积，得到的图像分别记为B⁺和I⁺；

第3.3步、根据图像B⁺和I⁺，从子区域中获取候选种子点；

第3.4步、从候选种子点中根据候选种子点间的距离剔除冗余的候选种子点，将剩余的候选种子点作为真实种子点。

5.根据权利要求3所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，第2.2步中，所述预设的范围为L_min～L_max，L_min=100±20像素，L_max=20000±200像素。

6.根据权利要求4所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，第3.3步中，对于子区域里面的每一个点o，如果该点的像素值B⁺(o)大于预设的阈值T3，而且I⁺(o)为局部极大值，将此点作为候选种子点，T3的取值范围为60～100个像素，其中，B⁺(o)表示点o的像素值，I⁺(o)表示点的局部极大值。

7.根据权利要求4所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，第3.4步的具体过程为，从候选种子点中选出I⁺(o)最大的那个点，将该点记为o₁，剔除与点o₁的欧拉距离小于预设的阈值T4的候选种子点，重复此操作直至候选种子点集为空集或重复至少10次，将剩余的候选种子点作为真实种子点，T4的取值为最小细胞的半径的0.8～0.9倍。

8.根据权利要求1所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，第4.1步中的式I转化为式II：

$\min {\left|\right|\underset{o\&Element;V}{\mathrm{\&Sigma;}}{B}_L\left(o\right)-\underset{o\&Element;V}{\mathrm{\&Sigma;}}\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}f(o,o_j,r_j)\left|\right|}^{1/3}+\mathrm{\&lambda;}{\mathrm{\&Sigma;}}_{j=1}^k{w_{j}r}_j$    式II

s.t.r_j≥0,j＝1,2,…,k

式中，w_j为权重，其中，res(o)为用所有球函数拟合B_L后没有被拟合到的部分，f(o,o_j,r_j)表示第j个球函数，其

表达式： $f(o,o_j,r_j)=\left\{\begin{array}{cc}1& \mathrm{if}|o-o_j|\&le;r_j\\ \exp (-\frac{{\left(\right|o-o_j|-r_j)}^2}{0.1})& \mathrm{otherwise}\end{array}\right.,$ o为BL上的像素点，

o_j为第j个种子点的坐标值，r_j为第j个种子点的半径。

9.根据权利要求8所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，利用梯度投影算法计算半径r_j，利用均值算法计算位置o_j。

10.根据权利要求1至9中任一所述的基于L1最小化模型的全自动细胞定位方法，其特征在于，所述细胞为神经元胞体。

Images