胞状物体通用分割算法Cellpose解析：开发篇
source link: https://qixinbo.info/2020/10/24/cellpose-2/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
上一篇介绍了如何安装和使用Cellpose，相当于将Cellpose当做一个开箱即用的软件。实际上Cellpose还是一个开源的代码库，开发者可以深入研究它的算法，并进行调用、修改和完善等。
本文尝试对Cellpose的运行机理做一个研究，包括它的标注数据的格式、神经网络的架构、神经网络的输入和输出等。
标注数据格式
假设有这么一张要分割的图像，大小为10像素乘以10像素（选择这么小的像素矩阵以方便打印和查看数值），背底为黑色，图像中间有一个白色圆盘，即：
它的像素矩阵为：
1
2
3
4
5
6
7
8
9
10
11
array([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 255., 255., 255., 0., 0., 0., 0.],
[ 0., 0., 255., 255., 255., 255., 255., 0., 0., 0.],
[ 0., 255., 255., 255., 255., 255., 255., 255., 0., 0.],
[ 0., 255., 255., 255., 255., 255., 255., 255., 0., 0.],
[ 0., 0., 255., 255., 255., 255., 255., 0., 0., 0.],
[ 0., 0., 0., 255., 255., 255., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],
dtype=float32)
对其进行手动标注，得到标注后的目标矩阵为：
向量场数据格式
Cellpose不是直接使用上面的标注数据作为神经网络的输入，而是要将这些标注转为向量场（这才是Cellpose最最重要的创新点，不是将物体的边缘轮廓直接作为目标拟合，而是将物体内部的场作为目标，非常巧妙），即热源在中心的热扩散场Flow Field，原理图如下：
该数据转换可以使用如下API函数：
1
2
3
4
5
import sys
sys.path.append('D:\\repos\\cellpose')
from cellpose import dynamics, io, plot
flow1 = dynamics.labels_to_flows([mask])
根据该函数的API介绍：
从这些数值可以看出，如果要实现可视化，还需要使用以下函数进行转换：
1
flow2 = plot.dx_to_circ([flow1[0][2], flow1[0][3]])
得到的就是RGB图像，即：
因为上面的样例是一张非常小的图，可以通过另外一张大的样例图来更直观地看一下效果：
向量场计算原理
上面概览了向量场的计算和可视化方式。这一节具体解析它的计算原理。
仍然以最开始的10乘10的小图为例，其mask为：
1
2
3
4
5
6
7
8
9
10
[[0 0 0 0 0 0 0 0 0 0]
[0 0 1 1 1 1 0 0 0 0]
[0 1 1 1 1 1 1 0 0 0]
[0 1 1 1 1 1 1 1 0 0]
[0 1 1 1 1 1 1 1 0 0]
[0 1 1 1 1 1 1 0 0 0]
[0 0 0 1 1 1 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]]
（1）计算向量场（cellpose中称为流场flow field），使用的是masks_to_flows这个函数：
首先寻找该mask中的物体：
1
slices = scipy.ndimage.find_objects(masks)
返回的是物体所在的行列切片：
1
slices = [(slice(1, 7, None), slice(1, 8, None))]
然后计算这些物体的中值直径大小：
1
dia = utils.diameters(masks)[0]
这里可以再插一句这个直径是怎样计算的：首先使用np.unique获得mask中特有元素的个数，元素个数就代表了每个物体的面积（unique函数去除其中重复的元素，并按元素由小到大返回一个新的无元素重复的元组或者列表），然后使用np.median计算这些面积的中值面积（如果有奇数个面积，则挑选中间的那个面积作为中值面积；如果有偶数个面积，则中间两个面积取平均作为中值面积），然后再求一个同样面积大小的圆的直径，即为中值直径。
接着计算在具体某个物体内它的质心的索引标识：
1
2
3
4
5
6
7
8
9
10
11
12
for i,si in enumerate(slices):
if si is not None:
sr,sc = si
ly, lx = sr.stop - sr.start + 1, sc.stop - sc.start + 1
y,x = np.nonzero(masks[sr, sc] == (i+1))
y = y.astype(np.int32) + 1
x = x.astype(np.int32) + 1
ymed = np.median(y)
xmed = np.median(x)
imin = np.argmin((x-xmed)**2 + (y-ymed)**2)
xmed = x[imin]
ymed = y[imin]
注意，这里的索引标识是以该物体内部自成一个坐标系，用到的关键一步就是代码中的nonzero那一行。然后同样使用np.median获取中位数，但因为不能确定该物体内是奇数还是偶数个像素（如果是偶数个像素，那么索引就是小数），所以先计算出中位数，再找出原来的索引序列中与该中位数最近的像素（即argmin那一行，通过距离判断），从而获得了真正的质心的索引标识。（这里还有一个伏笔：x和y在原有基础上又各加1，是为了后面计算扩散场时不从边界上开始计算）
接下来计算每个物体内的像素离其质心的距离，这种距离不是绝对距离，而是用之前的中位直径进行了某种归一化：
1
2
3
s2 = (.15 * dia)**2
d2 = (x-xmed)**2 + (y-ymed)**2
mu_c[sr.start+y-1, sc.start+x-1] = np.exp(-d2/s2)
这里的mu_c就是与原mask同样大小的距离变换图，因此在计算归一化距离时，虽然是相对于物体自己的质心的距离，但也通过sr.start的索引转换成了全局索引，即mu_c就是一幅以物体各自质心为中心的距离的一个个孤岛。
接下来计算热扩散场：
1
2
3
4
5
6
7
8
9
10
11
T = np.zeros((ly+2)*(lx+2), np.float64)
T = _extend_centers(T, y, x, ymed, xmed, np.int32(lx), niter)
def _extend_centers(T,y,x,ymed,xmed,Lx, niter):
for t in range(niter):
T[ymed*Lx + xmed] += 1
T[y*Lx + x] = 1/9. * (T[y*Lx + x] + T[(y-1)*Lx + x] + T[(y+1)*Lx + x] +
T[y*Lx + x-1] + T[y*Lx + x+1] +
T[(y-1)*Lx + x-1] + T[(y-1)*Lx + x+1] +
T[(y+1)*Lx + x-1] + T[(y+1)*Lx + x+1])
return T
T是一个ly+2乘以lx+2大小的一维数组，即它将原来的二维的物体区域给拉直成一个一维长条进行计算。
这里模拟的动作就是在质心处每次循环都投入一个数值为1的新的热源，然后计算整个物体内在这一热源下的热扩散，用的是八邻域计算。温度场在每次循环中都累加。
迭代循环结束后，还要叠加一个log型的温度场，这一步没有弄明白什么意思：
1
T[(y+1)*lx + x+1] = np.log(1.+T[(y+1)*lx + x+1])
接下来计算物体内的温度梯度场（即某点的x方向梯度就是左右两个邻居的温度的差，y方向梯度就是上下两个邻居的温度的差）：
1
2
dy = T[(y+1)*lx + x] - T[(y-1)*lx + x]
dx = T[y*lx + x+1] - T[y*lx + x-1]
并将其存入场变量mu中：
1
2
mu = np.zeros((2, Ly, Lx), np.float64)
mu[:, sr.start+y-1, sc.start+x-1] = np.stack((dy,dx))
注意mu的shape，dy和dx是分别存储的，分别存入mu的第0维和第1维。
最后将dx和dy进行标准化：
1
mu /= (1e-20 + (mu**2).sum(axis=0)**0.5)
即：先求(dy)^2和(dx)^2，然后通过sum求((dy)^2+(dx)^2)，然后对其开平方，最后将原值除以该值，得到：
1
2
3
dy/[((dy)^2+(dx)^2)**0.5]
和
dx/[((dy)^2+(dx)^2)**0.5]
（2）将流场与物体概率结合起来：
上面计算了流场后，在labels_to_flows中将该流场与是否是物体的概率结合起来：
1
2
3
veci = [masks_to_flows(labels[n][0])[0] for n in trange(nimg)]
flows = [np.concatenate((labels[n][[0]], labels[n][[0]]>0.5, veci[n]), axis=0).astype(np.float32)
for n in range(nimg)]
即这里的flows是四个量的结合：flows的形状是(4, Ly, Lx)，其中flows[k][0] 是标签labels[k]，flows[k][1]判断是否是细胞，所以值也就只有两个，0和1，其中在计算时判断是否是物体的阈值设为了0.5，flows[k][2]是Y方向的流场，flows[k][3]是X方向的流场。
从流场复原掩膜
上面解析了如何通过掩膜计算向量场（即流场）。流场是作为该算法中的输出，因此，算法推理后得到的输出也是流场形式，那么进一步地，需要根据流场反推出掩膜，以标识物体的形貌。
这里的分析要用到上面计算出来的流场flows。
为了更具通用性，在这一节中使用的掩膜中有两个物体，从而探究更一般的掩膜复原过程。
原始掩膜的矩阵为：
1
2
3
4
5
6
7
8
9
10
[[0 0 1 1 0 0 0 0 0 0]
[0 0 1 1 0 0 0 0 0 0]
[0 0 1 1 0 0 0 0 0 0]
[0 0 0 0 0 0 2 2 2 2]
[0 0 0 0 0 0 2 2 2 2]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]]
下面的步骤默认是不知道这个原始掩膜的，已知条件只有它的流场形式，通过上节的计算方法，得到该流场为（这里只显示第3和第4个元素，即X和Y方向的流场）：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
[[[ 0. 0. 0.7728456 0.7343627 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. -0.07701479 0.
0. 0. 0. 0. 0. ]
[ 0. 0. -0.8023549 -0.7343627 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0.6282049 0.9970669 0.72136474 0.26260668]
[ 0. 0. 0. 0. 0.
0. -0.67413443 -0.99975467 -0.80442435 -0.29906148]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]]
[[ 0. 0. 0.6345941 -0.67875725 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 1. -0.99702996 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0.59684724 -0.67875725 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0.778048 0.07653494 -0.69255537 -0.96490294]
[ 0. 0. 0. 0. 0.
0. 0.73860866 -0.0221485 -0.5940551 -0.9542338 ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. ]]]
接下来看怎样通过流场复原掩膜的。
（1）动力学计算
这一步在follow_flows函数中进行。该函数用到的参数只有一个与dP相关的数（还有一个迭代步数，有默认值）。
dP是流场flows的一部分，包含Y方向的流场和X方向的流场，可由下面的方式获得：
1
dP = flows[0][2:]
follow_flows函数所传入的实际参数在Cellpose中为：
1
-1 * dP * (cellprob > cellprob_threshold) / 5.
这其中有两个系数非常重要，一个是-1，另一个是5。
它们两者的作用分别是：-1是为了使符号反转，即原来正向热流，现在要使它往回流，即原流场是热源中的热往四周散，现在要复原了，所以要把热流反向，让热流都汇聚到热源上；5是为了将热流梯度值变小，这样能缓缓地流，保证热流流向临近点，防止一下冲到别的点上。
将这个表达式记为新dP，其值变为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
[[[-0. -0. -0.15456912 -0.14687255 -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. 0.01540296 -0.
-0. -0. -0. -0. -0. ]
[-0. -0. 0.16047098 0.14687255 -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0.12564097 -0.19941339 -0.14427295 -0.05252134]
[-0. -0. -0. -0. -0.
-0. 0.13482688 0.19995093 0.16088487 0.0598123 ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]]
[[-0. -0. -0.12691882 0.13575146 -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0.2 0.199406 -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0.11936945 0.13575146 -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0.1556096 -0.01530699 0.13851108 0.19298059]
[-0. -0. -0. -0. -0.
-0. -0.14772174 0.0044297 0.11881103 0.19084677]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. ]]]
首先根据dP的尺寸来生成动力学计算所在的网格（Cellpose能进行三维分割，不过这里只分析二维情形，即生成二维网格），即：
1
2
3
shape = np.array(dP.shape[1:]).astype(np.int32)
p = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]), indexing='ij')
p = np.array(p).astype(np.float32)
因为原图是一个10乘10的图像，那么这里的网格节点就是：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
[[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
[3. 3. 3. 3. 3. 3. 3. 3. 3. 3.]
[4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
[5. 5. 5. 5. 5. 5. 5. 5. 5. 5.]
[6. 6. 6. 6. 6. 6. 6. 6. 6. 6.]
[7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]
[8. 8. 8. 8. 8. 8. 8. 8. 8. 8.]
[9. 9. 9. 9. 9. 9. 9. 9. 9. 9.]]
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]]]
然后只计算流场梯度不为0的节点以加快计算速度，因此先得到这些不为0节点的索引：
1
inds = np.array(np.nonzero(np.abs(dP[0])>1e-3)).astype(np.int32).T
这些索引在该例中为：
1
2
3
4
5
6
7
8
9
10
11
12
13
[[0 2]
[0 3]
[1 3]
[2 2]
[2 3]
[3 6]
[3 7]
[3 8]
[3 9]
[4 6]
[4 7]
[4 8]
[4 9]]
然后执行以下函数：
1
p = steps2D(p, dP, inds, niter)
具体它干的工作就是不断地在格点上累加热流。
结果就是：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
[[[0. 0. 1.0098284 1.1413077 0. 0. 0.
0. 0. 0. ]
[1. 1. 1. 1.0362048 1. 1. 1.
1. 1. 1. ]
[2. 2. 1.1182916 1.1164871 2. 2. 2.
2. 2. 2. ]
[3. 3. 3. 3. 3. 3. 4.105441
3.8231199 4.041441 4.0182114]
[4. 4. 4. 4. 4. 4. 3.8856084
3.9896855 3.9375546 3.935263 ]
[5. 5. 5. 5. 5. 5. 5.
5. 5. 5. ]
[6. 6. 6. 6. 6. 6. 6.
6. 6. 6. ]
[7. 7. 7. 7. 7. 7. 7.
7. 7. 7. ]
[8. 8. 8. 8. 8. 8. 8.
8. 8. 8. ]
[9. 9. 9. 9. 9. 9. 9.
9. 9. 9. ]]
[[0. 1. 2.8871436 3.0260515 4. 5. 6.
7. 8. 9. ]
[0. 1. 2. 3.1274133 4. 5. 6.
7. 8. 9. ]
[0. 1. 2.811595 3.157068 4. 5. 6.
7. 8. 9. ]
[0. 1. 2. 3. 4. 5. 7.9200172
7.898424 7.9439383 7.9862995]
[0. 1. 2. 3. 4. 5. 7.990998
7.884873 7.924165 7.893717 ]
[0. 1. 2. 3. 4. 5. 6.
7. 8. 9. ]
[0. 1. 2. 3. 4. 5. 6.
7. 8. 9. ]
[0. 1. 2. 3. 4. 5. 6.
7. 8. 9. ]
[0. 1. 2. 3. 4. 5. 6.
7. 8. 9. ]
[0. 1. 2. 3. 4. 5. 6.
7. 8. 9. ]]]
p的值域也是有规则的，即大于0且小于原图尺寸，因为：
1
p[0,y,x] = min(shape[0]-1, max(0, p[0,y,x] - dP[0,p0,p1]))
单看p中具体数值，没有意义，应该看p现在的值与之前的网格格点标识之间的差值，这才是有意义的，表明在该点上的热流累加（或散失），由此来标识经过动力学计算后热源所在位置。
另外剧透一下，从下面的分析可以看出，目前p中的值经过整型化后，与原先的格点索引不同的位置上的数值正是热源所在位置，比如第一排第三个元素的原先索引是[0, 2]，现在变成了[1, 2]，那么这个[1, 2]正是热源所在位置。
（2）复原掩膜：
1
2
3
for i in range(dims):
pflows.append(p[i].flatten().astype('int32'))
edges.append(np.arange(-.5-rpad, shape0[i]+.5+rpad, 1))
pflows是对浮点数的p展平后进行整型化，即：
1
2
3
4
5
6
7
8
9
10
[array([0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,
1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 3, 4, 4, 4, 4, 4, 4,
4, 4, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9]),
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1,
2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 7, 7, 7, 7, 0, 1, 2, 3,
4, 5, 7, 7, 7, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5,
6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7,
8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9])]
而edges是直方图中的各个区间范围，其中的参数rpad控制了直方图的边缘填充大小，比如rpad如果为0，那么直方图的大小就是与原图大小相等，即周围没有填充0，而如果设置rpad为非零，比如设rpad为1，那么：
1
2
3
4
[array([-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5,
9.5, 10.5]),
array([-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5,
9.5, 10.5])]
而因为pflows中的值都是大于0且小于原图尺寸的（step2D中的操作），即-1.5和-0.5之间以及9.5和10.5之间是没有元素的，就相等于在直方图周围填充了一圈0的边界，如果rpad为更多，比如默认为20，那么填充的边界宽度则为20。
这里将rpad设为0，即不要填充边界，那么：
1
2
edges = [array([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5]),
array([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5])]
这样以0.5为后缀，是为了方便取得下面的整型，比如0.5到1.5之间就是取的1。
然后就是至关重要的直方图计算环节：
1
h,_ = np.histogramdd(tuple(pflows), bins=edges)
这一步就是根据edges区间统计这些区间内相应流场强度，先从结果来向前推比较好理解，这里h.T（注意是将原来的h进行了以下转置，方便描述）的值为：
1
2
3
4
5
6
7
8
9
10
[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[0. 3. 0. 1. 1. 1. 1. 1. 1. 1.]
[0. 3. 0. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 0. 0. 1. 1. 1. 1. 1.]
[1. 1. 1. 5. 3. 1. 1. 1. 1. 1.]
[1. 1. 1. 0. 0. 1. 1. 1. 1. 1.]
[1. 1. 1. 0. 0. 1. 1. 1. 1. 1.]]
第一个元素就代表当y方向的pflows为-0.5到0.5之间（就是0）时x方向的pflows为-0.5到0.5之间（就是0）的数的个数，那么看pflows就是(0, 0)这个元素；第二个元素就代表y方向的pflows为0时x方向的pflows为1的数的个数，发现就是(1, 0)这一对。比如第三行的3代表当y方向的pflows为1.5到2.5之间（就是2）时，x方向的pflows为0.5到1.5之间（就是1）的数的个数是3个，分别在第3个、第13个、第23个元素对。
这部分的理解可以根据下面更通俗的一个例子来辅助。
np.histgram就是在一维上看数据落在bins上的个数。
np.histgram2d就是在二维上看数据落在bins上的个数，以如下例子为例：
1
2
3
4
5
6
7
8
import numpy as np
xedges = [0, 1, 3, 5]
yedges = [0, 2, 3, 4, 6]
x = np.random.normal(2, 1, 10)
y = np.arange(10)
H, xedges, yedges = np.histogram2d(x, y, bins=(xedges, yedges))
上述代码中xedges定义了x维度上的值所划分的范围（即统计0到1之间（注意不含1）的值的个数、1到3之间（注意不含3）的值的个数、3到5之间（注意这里包含5）的值的个数），yedges定义了y维度上的值所划分的范围，x是一个含有10个以2为期望、以1为标准差的正态分布的数值：
1
2
array([1.4965097 , 0.46480962, 1.65365048, 1.56658642, 2.05414053,
3.209907 , 2.87861765, 0.17336843, 2.97341494, 2.12467305])
y就是一个从0到9依次递增的10个数值：
1
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
而得到的H为：
1
2
3
array([[1., 0., 0., 0.],
[1., 1., 1., 2.],
[0., 0., 0., 1.]])
将其转置一下：
1
2
3
4
5
>>> H.T
array([[1., 1., 0.],
[0., 1., 0.],
[0., 1., 0.],
[0., 2., 1.]])
H.T是一个4乘3的矩阵，第一个元素代表x维度上0到1之间（注意不含1）且y维度上0到2之间（注意不含2）的数的个数，即(0.46480962, 1)这个数，第二个元素代表x维度上1到3之间（注意不含3）且y维度上0到2之间（注意不含2）的数的个数，即(1.4965097, 0)这个数，第三个元素代表x维度上3到5之间（注意!!这里包含5）且y维度上0到2之间（注意不含2）的数的个数，没有这样的数，所以为0。看一下H.T中倒数第二个元素，它代表x维度上1到3之间（注意不含3）且y维度上4到6之间（注意!!包含6）的数的个数，即(2.05414053, 4)和(2.87861765, 6)这两个数。因此，在edges范围两端的值，都是包含在范围内，而在内部的值，左边值包含，而右边值不包含。
对于bins这个参数，上面是指定的范围序列，也可以指定一个整数值，表示分成多少份。
以上是关于那个直方图矩阵的理解，对于其实际物理意义，可以想象成这样：如果pflows还是保持原来的格点索引，即：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
[[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
[3. 3. 3. 3. 3. 3. 3. 3. 3. 3.]
[4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
[5. 5. 5. 5. 5. 5. 5. 5. 5. 5.]
[6. 6. 6. 6. 6. 6. 6. 6. 6. 6.]
[7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]
[8. 8. 8. 8. 8. 8. 8. 8. 8. 8.]
[9. 9. 9. 9. 9. 9. 9. 9. 9. 9.]]
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]]]
那么上面的统计矩阵上的元素应该都为1，可以想象成这个棋盘上每个格点上都有一粒米，那么再看pflows经过了动力学计算后的形式：
1
2
3
4
5
6
7
8
9
10
[array([0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,
1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 3, 4, 4, 4, 4, 4, 4,
4, 4, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9]),
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1,
2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 7, 7, 7, 7, 0, 1, 2, 3,
4, 5, 7, 7, 7, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5,
6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7,
8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9])]
以及再贴一遍直方图统计结果：
1
2
3
4
5
6
7
8
9
10
[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[0. 3. 0. 1. 1. 1. 1. 1. 1. 1.]
[0. 3. 0. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 0. 0. 1. 1. 1. 1. 1.]
[1. 1. 1. 5. 3. 1. 1. 1. 1. 1.]
[1. 1. 1. 0. 0. 1. 1. 1. 1. 1.]
[1. 1. 1. 0. 0. 1. 1. 1. 1. 1.]]
看pflows，比如第一行，x维度上两个0变成了两个1，表示x维度上原来0号标识的位置上的米移到了1号标识上，y维度保持不变，对应于直方图统计就是将两粒米右移一格，x维度上有两个2也变成了1，意味着这两个x维度上2号标识上的米左移到了1号标识上。同理，x、y坐标为(3, 6)的点现在变成了(4, 7)，也是把这个米从这两个坐标点上进行了移动。
上面是以相对位置角度来理解，另一个角度来看，从绝对位置来看，pflows中的要移动米粒的位置上面的值正是热源所在坐标位置：还是看x维度上的那两个变为1标识的0标识，它原来的标识是[0, 2]，现在变成了[1, 2]，这[1, 2]正是热源位置（下面的求解可知道）。
关于直方图统计这一部分理解结束！
下面接着看。
得到直方图统计结果后，对其进行一个核为5的最大化池化（池化的目的是为了滤波，去掉噪声）：
1
2
for i in range(dims):
hmax = maximum_filter1d(hmax, 5, axis=i)
1
2
3
4
5
6
7
8
9
[[3. 3. 3. 3. 3. 3. 1. 1. 1. 1.]
[3. 3. 3. 3. 3. 5. 5. 5. 5. 5.]
[3. 3. 3. 3. 3. 5. 5. 5. 5. 5.]
[3. 3. 3. 3. 3. 5. 5. 5. 5. 5.]
[1. 1. 1. 1. 1. 5. 5. 5. 5. 5.]
[1. 1. 1. 1. 1. 5. 5. 5. 5. 5.]
[1. 1. 1. 1. 1. 3. 3. 3. 3. 3.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
然后找出原直方图统计中局部最大值的位置：
1
seeds = np.nonzero(np.logical_and(h-hmax>-1e-6, h>2))
具体判断标准有两个：一是看原值与经过池化后的值的差值，如果差值为0，则表明这个位置是局部极大值；二是该局部极大值需要大于2（注意，原代码中这个地方是10，这里因为原图只有10乘10，所以改了一下，变为2，这样才能检测出来，所以这个10也是一个可调参数），即它的热流汇聚要足够大，是一个大的热源，也是为了排除小的干扰。
seeds的值为（即局部极大值的坐标）：
1
seeds = (array([1, 1, 3], dtype=int64), array([2, 3, 7], dtype=int64))
seeds的第0个元素是局部极大值的x坐标，第1个元素是其对应的y坐标。下面将这些坐标成对提取出来：
1
pix = list(np.array(seeds).T)
1
pix = [array([1, 2], dtype=int64), array([1, 3], dtype=int64), array([3, 7], dtype=int64)]
然后构建一个3乘3（注意这是二维，三维下是3乘3乘3）的窗口：
1
expand = np.nonzero(np.ones((3,3)))
这个窗口内的数值是节点索引：
1
expand = (array([0, 0, 0, 1, 1, 1, 2, 2, 2], dtype=int64), array([0, 1, 2, 0, 1, 2, 0, 1, 2], dtype=int64))
接下来是两层循环。先看内部循环，它是对pix的长度进行循环（再次明确一下，pix是局部极大值的坐标索引）：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
for k in range(len(pix)):
if iter==0:
pix[k] = list(pix[k])
newpix = []
iin = []
for i,e in enumerate(expand):
epix = e[:,np.newaxis] + np.expand_dims(pix[k][i], 0) - 1
epix = epix.flatten()
iin.append(np.logical_and(epix>=0, epix<shape[i]))
newpix.append(epix)
iin = np.all(tuple(iin), axis=0)
newpix = tuple(newpix)
igood = h[newpix]>2
for i in range(dims):
pix[k][i] = newpix[i][igood]
if iter==4:
pix[k] = tuple(pix[k])
这里面的变量较多，各自具体的意义为：pix刚开始是局部极大值的坐标索引，epix是在某一个维度上expand滑动窗加上pix索引减去1，newpix是在两个维度上对epix进行合成，这样来说，newpix就是以pix中的坐标为中心的九个像素坐标索引，比如pix中某个像素是[1, 2]，那么newpix就是：
1
newpix = [array([0, 0, 0, 1, 1, 1, 2, 2, 2], dtype=int64), array([1, 2, 3, 1, 2, 3, 1, 2, 3], dtype=int64)]
然后看之前的直方图统计矩阵h在这九个像素点上是否大于2，即为igood的值，即是否有热流流入，以上面的[1, 2]中心像素所产生的九个像素为例，igood就是：
1
igood = [False False False False True True False False False]
即仍然是[1, 2]和[1, 3]中的热流强度大于2，然后把这两个坐标提取出来再次放入pix中，即下面这句：
1
2
for i in range(dims):
pix[k][i] = newpix[i][igood]
此时pix变为：
1
pix = [[array([1, 1], dtype=int64), array([2, 3], dtype=int64)], array([1, 3], dtype=int64), array([3, 7], dtype=int64)]
对比一下pix最开始的坐标索引，可以发现pix中第一个元素由原来的局部极大值索引[1, 2]，变成了以它为中心的九个像素中直方图统计大于2的像素坐标索引。
那么对原来pix中所有极大值索引都循环一遍后，得到了新的pix：
1
2
3
[[array([1, 1], dtype=int64), array([2, 3], dtype=int64)],
[array([1, 1], dtype=int64), array([2, 3], dtype=int64)],
[array([3, 4], dtype=int64), array([7, 7], dtype=int64)]]
以上就是内层循环的作用：以最开始的局部极大值所在像素为中心，在其上罩一个3乘3的窗口，看该窗口内直方图统计是否大于2，如果是，就将该邻居像素的坐标索引加入pix中。
那么外层循环就是不断地重复这个内层循环5次（这里的5也是一个可调参数），即一步步地查找是否由邻居像素的热流大于2。
内层循环和外层循环加起来实现的效果就是：模仿漫水填充，查找最开始的最大热源所辐射的范围，看哪些邻居像素属于该热源。
在该例中，经过上述内外循环后，pix的值为：
1
2
3
4
5
6
7
8
9
10
11
12
[(array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int64),
array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3], dtype=int64)),
(array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int64),
array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3], dtype=int64)),
(array([3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4], dtype=int64),
array([7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
7, 7, 7, 7, 7, 7, 7, 7, 7, 7], dtype=int64))]
虽然看起来里面有很多元素，但pix的length只有3，表明这些元素都是由最开始的3个局部极大值所辐射的；另一方面，很多元素都是重复的，实际只有[1, 2]、[1, 3]、[3, 7]、[4, 7]这四个元素。
然后对这些元素进行label：
1
2
3
M = np.zeros(h.shape, np.int32)
for k in range(len(pix)):
M[pix[k]] = 1+k
可得M矩阵的数值（注意它这里是h.shape，即M是与直方图统计所对应的，不是原图）：
1
2
3
4
5
6
7
8
9
10
[[0 0 0 0 0 0 0 0 0 0]
[0 0 2 2 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 3 0 0]
[0 0 0 0 0 0 0 3 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]]
从上面已经知道，pflows中异常值就代表了热源位置，所以这里根据pflows的索引，得到其上的label值（即将pflows传入M）：
1
2
3
4
for i in range(dims):
pflows[i] = pflows[i] + rpad
M0 = M[tuple(pflows)]
所以，与原图对应的label值（就是掩膜）就是：
1
2
3
[0 0 2 2 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0 0 0 3
3 3 3 0 0 0 0 0 0 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
这里得到掩膜后，中间还插了一段去掉非常大的掩膜的代码：
1
2
3
4
_,counts = np.unique(M0, return_counts=True)
big = np.prod(shape0) * 0.4
for i in np.nonzero(counts > big)[0]:
M0[M0==i] = 0
因为这个例子的掩膜比较小，因此这段代码没起作用。
从目前掩膜M0的值来看，它里面的数值是2和3，没有从1开始。
下面就是转换一下（这里用到的技术就是np.unique中的return_inverse这个参数，是返回这些unique数的索引，非常巧妙）：
1
2
_,M0 = np.unique(M0, return_inverse=True)
M0 = np.reshape(M0, shape0)
此时M0变为：
1
2
3
[0 0 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 2
2 2 2 0 0 0 0 0 0 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
最后再变成原图的尺寸即可：
1
M0 = np.reshape(M0, shape0)
1
2
3
4
5
6
7
8
9
10
[[0 0 1 1 0 0 0 0 0 0]
[0 0 1 1 0 0 0 0 0 0]
[0 0 1 1 0 0 0 0 0 0]
[0 0 0 0 0 0 2 2 2 2]
[0 0 0 0 0 0 2 2 2 2]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]]
UnetModel模型类
Cellpose的算法是基于Unet架构的。
基础模型调用在：
1
2
3
4
5
6
7
class UnetModel():
self.net = resnet_style.CPnet(nbase, nout=self.nclasses,
residual_on=residual_on,
style_on=style_on,
concatenation=concatenation)
self.net.hybridize(static_alloc=True, static_shape=True)
self.net.initialize(ctx = self.device)
即，UnetModel这个类实际是CPnet类的一个实例化对象，因此，CPnet类才是最基础的算法架构所在：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
class CPnet(gluon.HybridBlock):
def __init__(self, nbase, nout, residual_on=True, style_on=True, concatenation=False, **kwargs):
super(CPnet, self).__init__(**kwargs)
with self.name_scope():
self.nbase = nbase
self.downsample = downsample(nbase, residual_on=residual_on)
self.upsample = upsample(nbase, residual_on=residual_on, concatenation=concatenation)
self.output = batchconv(nout, 1)
self.make_style = make_style()
self.style_on = style_on
def hybrid_forward(self, F, data):
#data = self.conv1(data)
T0 = self.downsample(data)
style = self.make_style(T0[-1])
style0 = style
if not self.style_on:
style = style * 0
T0 = self.upsample(style, T0)
T0 = self.output(T0)
return T0, style0
CPnet基于MXNet/Gluon（Gluon是MXNet的动态图接口），这里还使用了MXNet的混合式编程方法，即继承自gluon.HybridBlock类。混合式编程是将命令式编程与符号式编程混合在一起，既可以通过命令式编程使代码编写起来更容易和直观，也可以通过符号式编程来获得更好的计算性能和可移植性。MXNet通过HybridSequential类和HybridBlock类构建的模型可以调用hybridize函数将命令式程序转成符号式程序。具体教程见下方链接：
命令式和符号式混合编程
CPnet类有两种算法结构可以选择，一种是经典的Unet结构，一种是以ResNet为backbone的Unet结构，两者可以通过residual_on这个参数来开关。
CellposeModel模型类
CellposeModel类是基于上面的UnetModel类，但两者也有不同：UnetModel是通用的算法架构，其训练集的data和label是通用的、未经特殊处理的；而CellposeModel的训练集的label实际是在原mask基础上计算的flow field，即定制化的label。
主要区别在这里：
对于UnetModel，其label为：
1
2
3
4
5
6
7
if self.nclasses==3:
print('computing boundary pixels')
train_classes = [np.stack((label, label>0, utils.distance_to_boundary(label)), axis=0).astype(np.float32)
for label in tqdm(train_labels)]
else:
train_classes = [np.stack((label, label>0), axis=0).astype(np.float32)
for label in tqdm(train_labels)]
对于CellposeModel，其label为：
1
train_flows = dynamics.labels_to_flows(train_labels, files=train_files)
SizeModel模型类
SizeModel类是用来估计图像中物体尺寸（中位直径）的模型。
Cellpose类
Cellpose类是对CellposeModel和SizeModel的整合，即对这两个模型的统一调用。
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK