胞状物体通用分割算法Cellpose解析:开发篇

概览

上一篇介绍了如何安装和使用Cellpose,相当于将Cellpose当做一个开箱即用的软件。实际上Cellpose还是一个开源的代码库,开发者可以深入研究它的算法,并进行调用、修改和完善等。
本文尝试对Cellpose的运行机理做一个研究,包括它的标注数据的格式、神经网络的架构、神经网络的输入和输出等。

标注数据格式

假设有这么一张要分割的图像,大小为10像素乘以10像素(选择这么小的像素矩阵以方便打印和查看数值),背底为黑色,图像中间有一个白色圆盘,即:
demo
它的像素矩阵为:

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)

对其进行手动标注,得到标注后的目标矩阵为:
mask
即,白色区域标为1,背景区域标为0(这里因为是手动标注,得到的标注不一定非常严格准确)。如果是标注了多个区域,则该mask矩阵里的标注值依次递增,即1、2、3等等。(注意这个地方一定要通过skimage.io.imread读取,不要用opencv,否则读取的是错误的值)

向量场数据格式

Cellpose不是直接使用上面的标注数据作为神经网络的输入,而是要将这些标注转为向量场(这才是Cellpose最最重要的创新点,不是将物体的边缘轮廓直接作为目标拟合,而是将物体内部的场作为目标,非常巧妙),即热源在中心的热扩散场Flow Field,原理图如下:

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介绍:
flow-api
可知,返回的flow1中包含了向量场的Y分量和X分量:
flow-Y
flow-X

从这些数值可以看出,如果要实现可视化,还需要使用以下函数进行转换:

1
flow2 = plot.dx_to_circ([flow1[0][2], flow1[0][3]])

得到的就是RGB图像,即:
flow-visual

因为上面的样例是一张非常小的图,可以通过另外一张大的样例图来更直观地看一下效果:
contrast

向量场计算原理

上面概览了向量场的计算和可视化方式。这一节具体解析它的计算原理。
仍然以最开始的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的整合,即对这两个模型的统一调用。