neuroGPS算法介绍

NeuroGPS（neuro globel position system）是由华中科技大学全廷伟教授于2013年发表于scientific reports用来实现神经元胞体定位的方法，该方法基于L1范数可实现对小鼠脑神经元的形状、大小、密度的广泛多样性具有很强鲁棒性，该方法能够消除复杂神经元突起，特别是粗大的树突对定位的干扰。
原论文：NeuroGPS automated localization of neurons for brain circuits using L1 minimization model
以下为本人对论文理解以及整理：

模型建立

首先，通过对原始图像进行二值化、腐蚀操作后，得到二值化图像BL，假设神经元胞体为一个理想的球体，也就是说对于胞体公式，可以用f(o,oi,ri)这个球方程表示，其中oi和ri分别表示第i个胞体坐标和半径，o表示体素V中坐标；对于体素V中非胞体数据，则统一用res（o）表示，所以BL即等于球体和残差部分总和，即：

而对于球方程f(o,oi,ri)，由下面方程确定：

对于球内的像素，直接为1，之外像素值大小用一个幂函数确定。
对于方程1，需要估计的未知数即球中心坐标和球半径大小，使得方程1拟合的最好，也即是尽量使得res（o）的大小最小，最经典的算法就是最小二乘法，即求下面式子最小值：

其中|| ||表示2范数平方，V是BL中体积像素的坐标集。
由于考虑到神经元的空间分布是稀疏的，即神经元胞体之间空间分布是不重叠的，如果我们为某个数据集设置了许多潜在的胞体位置，这些位置中的大多数应该是假的，并且相关半径应该是零，也就是半径具有稀疏性，这个可以利用L1范数规则化公式3，提高半径稀疏性，即：
其中λ是球函数的参数拟合误差与球函数半径的稀疏性之间的折衷，为了保持相同的物理维度，目标函数的第一项被修改为误差的2范数平方的1/3次方。
后面就是对上述方程4的求解问题，作者分别使用梯度投影法和平均法更新半径和位置。最后对于所有最优半径，如果大于给定的阈值，则将相应的位置视为有效位置。

具体实现步骤

1、图像二值化处理

图像二值化处理的作用其实就是考虑到神经元的稀疏分布，如果对整个图像进行分析，太过复杂，而如果只对神经元那部分区域分析则将显著提高计算效率二值化操作就是为了后续子区域的提取（我觉得这是一个创新点），二值化的图像B采用下面这个公式：

其中I(o)和C(o)分别表示原始图像和背景图像中的体像素坐标o的灰度值，背景图像C（o）通过下面步骤获取：
步骤一：计算i*（o）=min（I(o)，thre），其中thre是用于区分前景和背景的粗略估计值；
步骤二：将i*与平均模板（9×9×1）卷积20次，得到背景图像C（o）。

2、侵蚀操作

这一步，就是子区域提取步骤，只不过这个步骤是通过侵蚀操作实现，具体如下步骤：
步骤一：侵蚀二值化图像B，考虑一个333区域内，如果该区域内所有像素值的和小于一个阈值T，就将该点值设为0；
步骤二：经过上一步后，提取并分析体积像素数为100~20000的连通区，提取完这部分后就将提取区域中的像素值设定为0,；
步骤三：从步骤1和步骤2一直重复更新B，直到B中任意连接的区域体像素小于100。
作者在第一次侵蚀操作中，侵蚀阈值T被设置为9，并且对于后续的侵蚀，侵蚀阈值T以0.027的步长持续增加，而对于提取区域的有效体积像素数目（上诉步骤中的100~20000）可以由提取区域中的最小神经元大小、体积像素大小和最大预定神经元数目来确定，并由NeuroGPS自动估计。

3、种子点选取

将之前提取到的子区域的二值化和原始信号分别嵌入到新的空长方体区域中，用BL和IL表示，然后将其都与为1的模板（775）卷积，生成对应的滤波信号B+和I+，将满足B+（o）大于70且I+ (o)为局部峰值的位置o作为候选种子。冗余种子点的删除，从候选种子中选择I+ (o)值最高的种子o1，删除满足位置o与位置o1的欧几里德距离小于4.8 µm的候选种子o。重复此过程，直到候选种子集为空或提取的种子数大于10。

4、方程4优化

得到候选种子点后，就是求解方程4的问题：
步骤1：初始化权重参数wi=1，i=1,2,…,k（这里是引用了“更新加权的L1最小算法”，可保证求解方程4时自动修改稀疏性的强度，在λ值较小的情况下也能提供强稀疏性）；
步骤2：解决优化问题：（注意方程式中BL，表示操作的是被提取出来的子区域二值化值）

（1）设置种子坐标和半径初始值（o1t，r1t），（o2t，r2t），…，（okt，rkt），t=0；（这里估计对于半径就是自行设定一个统一初始值）
（2）使用梯度投影法更新半径r；
（3）使用平均法更新位置：

其中，就是在第i个这个球内部像素点原始灰度大小，（注意：这里球表示为ot和rt+1，不是rt；还有这里用的是Ii，表示用的是灰度值）；
一直重复（2）、（3）步骤，直到满足终止条件；
步骤3：更新权重：

步骤4：一直重复步骤2、3，直到所有参数收敛（提示：步骤2是为了获得收敛的o和r，步骤3是为了获得收敛的w值）。

5、有效位置合并

经过前面步骤得到种子点的新的半径值，通过设定一个阈值Lmin剔除掉半径小于这个阈值的球，然后如果两个有效位置之间的距离小于这两个位置半径之和的70%，则将这两个位置合并为一个位置，否则保持不变。最后合并或解合后的位置作为识别出的神经元位置。

结语

以上为neuroGPS算法全部过程，如果想实现应该不难，等我python实现后再更新，不足之处请批评指正。