概览
上一篇介绍了如何安装和使用Cellpose，相当于将Cellpose当做一个开箱即用的软件。实际上Cellpose还是一个开源的代码库，开发者可以深入研究它的算法，并进行调用、修改和完善等。
本文尝试对Cellpose的运行机理做一个研究，包括它的标注数据的格式、神经网络的架构、神经网络的输入和输出等。
标注数据格式
假设有这么一张要分割的图像，大小为10像素乘以10像素（选择这么小的像素矩阵以方便打印和查看数值），背底为黑色，图像中间有一个白色圆盘，即：
它的像素矩阵为：1
2
3
4
5
6
7
8
9
10
11array([[ 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)
对其进行手动标注，得到标注后的目标矩阵为：
即，白色区域标为1，背景区域标为0（这里因为是手动标注，得到的标注不一定非常严格准确）。如果是标注了多个区域，则该mask矩阵里的标注值依次递增，即1、2、3等等。（注意这个地方一定要通过skimage.io.imread读取，不要用opencv，否则读取的是错误的值）
向量场数据格式
Cellpose不是直接使用上面的标注数据作为神经网络的输入，而是要将这些标注转为向量场（这才是Cellpose最最重要的创新点，不是将物体的边缘轮廓直接作为目标拟合，而是将物体内部的场作为目标，非常巧妙），即热源在中心的热扩散场Flow Field，原理图如下：
该数据转换可以使用如下API函数：1
2
3
4
5import sys
sys.path.append('D:\\repos\\cellpose')
from cellpose import dynamics, io, plot
flow1 = dynamics.labels_to_flows([mask])
根据该函数的API介绍：
可知，返回的flow1中包含了向量场的Y分量和X分量（这里的Y和X是相对于图像坐标系来说，即垂直为Y，水平为X，但实际上在numpy计算时，约定是垂直为x，水平为y，所以两者是一个颠倒的关系，注意分辨）：
从这些数值可以看出，如果要实现可视化，还需要使用以下函数进行转换：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
12for 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
3s2 = (.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
2dy = T[(y+1)*lx + x] - T[(y-1)*lx + x]
dx = T[y*lx + x+1] - T[y*lx + x-1]
并将其存入场变量mu中：1
2mu = 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
3dy/[((dy)^2+(dx)^2)**0.5]
和
dx/[((dy)^2+(dx)^2)**0.5]
（2）将流场与物体概率结合起来：
上面计算了流场后，在labels_to_flows中将该流场与是否是物体的概率结合起来：1
2
3veci = [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个元素，即图像坐标系的Y和X方向的流场）：
1 | [[[ 0. 0. 0.7728456 0.7343627 0. |
通过流场复原掩膜的函数如下：1
mask, p = dynamics.compute_masks(flow[0][2:], flow[0][1], flow_threshold=0.0, min_size=1, interp=False)
得到的掩膜为（因为本例的图像非常小，需要设置如上的参数，及在代码中计算seeds时还要改一处硬编码的参数，具体见下面解析）：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]]
可以发现能准确地复原产生该流场的掩膜。
下面具体看一下怎样计算的。
根据可能性筛选有效像素
（0）根据probability计算有效像素1
cp_mask = cellprob > cellprob_threshold
这一步是通过probability与其阈值（默认是0）进行比较，筛选出有效的像素。
重整流场强度
接下来进入follow_flows函数：1
2p, inds = follow_flows(dP * cp_mask / 5., niter=niter, interp=interp,
use_gpu=use_gpu, device=device)
该函数用到的重要参数只有一个与dP相关的数（还有一个迭代步数，有默认值）。
dP是流场flows的一部分，包含Y方向的流场和X方向的流场，可由下面的方式获得：1
dP = flows[0][2:]
follow_flows函数所传入的实际参数为：1
dP * cp_mask / 5.
即将原来的流场强度经过probability的选择后，再除以5，即降低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
3shape = 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.astype(np.float32), inds, niter)
具体它干的工作就是不断地在网格点阵上累加热流，即进行欧拉积分（步长为1）。
具体代码解释如下：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
34def steps2D(p, dP, inds, niter):
# 获得图像的尺寸，该例中就是(10, 10)
shape = p.shape[1:]
print("dp = ", dP)
# 进行niter次迭代积分
for t in range(niter):
# 只对有效像素进行循环，即只看inds索引，这样提高计算速度
for j in range(inds.shape[0]):
# 取出有效索引inds中的y和x坐标
y = inds[j,0]
x = inds[j,1]
print("y = ", y)
print("x = ", x)
# 将上述坐标映射回网格点阵中的坐标
# 这里不能直接使用y和x
# 而是要经过一个整型化，即选择该像素最近的、且是热量来源方向的网格格点
p0, p1 = int(p[0,y,x]), int(p[1,y,x])
print("p0 = ", p0)
print("p1 = ", p1)
# 根据网格点阵的坐标，取得该像素点的y和x方向的梯度
step = dP[:,p0,p1]
print("step = ", step)
# 对图像的两个维度进行循环
for k in range(p.shape[0]):
print("p[k,y,x] 000000 = ", p[k,y,x])
# 使用欧拉积分：y_(t+1) = y_t + delta_t*dP
# 这里的步长delta_t设置为1，所以y_(t+1) = y_t +dP
p[k,y,x] = min(shape[k]-1, max(0, p[k,y,x] + step[k]))
print("p[k,y,x] 111111 = ", p[k,y,x])
print("p = ", p)
print("********")
return p
我们先看迭代一次的结果，此时网格点阵p变成了：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. ]
[1. 1. 1. 0.984597 1. 1.
1. 1. 1. 1. ]
[2. 2. 1.839529 1.8531275 2. 2.
2. 2. 2. 2. ]
[3. 3. 3. 3. 3. 3.
3.1256409 3.1994133 3.144273 3.0525212 ]
[4. 4. 4. 4. 4. 4.
3.865173 3.800049 3.8391151 3.9401877 ]
[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.1269188 2.8642485 4. 5.
6. 7. 8. 9. ]
[0. 1. 2. 2.800594 4. 5.
6. 7. 8. 9. ]
[0. 1. 2.1193695 2.8642485 4. 5.
6. 7. 8. 9. ]
[0. 1. 2. 3. 4. 5.
6.1556096 7.015307 7.861489 8.807019 ]
[0. 1. 2. 3. 4. 5.
6.147722 6.99557 7.881189 8.809154 ]
[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. ]]]
即原先的网格点阵索引都累加了该点位置上的流场梯度，即$y_{t+1}=y_t+dP$（步长取为1)。
再经过一次迭代，此时仍然是累加热流，但与上次稍有不同，看以下具体的中间的输出：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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158y = 0
x = 2
p0 = 0
p1 = 2
step = [0.15456912 0.12691882]
p[k,y,x] 000000 = 0.1545691192150116
p[k,y,x] 111111 = 0.3091382384300232
p[k,y,x] 000000 = 2.1269187927246094
p[k,y,x] 111111 = 2.2538375854492188
y = 0
x = 3
p0 = 0
p1 = 2
step = [0.15456912 0.12691882]
p[k,y,x] 000000 = 0.14687255024909973
p[k,y,x] 111111 = 0.30144166946411133
p[k,y,x] 000000 = 2.864248514175415
p[k,y,x] 111111 = 2.9911673069000244
y = 1
x = 3
p0 = 0
p1 = 2
step = [0.15456912 0.12691882]
p[k,y,x] 000000 = 0.9845970273017883
p[k,y,x] 111111 = 1.1391661167144775
p[k,y,x] 000000 = 2.8005940914154053
p[k,y,x] 111111 = 2.9275128841400146
y = 2
x = 2
p0 = 1
p1 = 2
step = [0. 0.2]
p[k,y,x] 000000 = 1.839529037475586
p[k,y,x] 111111 = 1.839529037475586
p[k,y,x] 000000 = 2.1193695068359375
p[k,y,x] 111111 = 2.3193695545196533
y = 2
x = 3
p0 = 1
p1 = 2
step = [0. 0.2]
p[k,y,x] 000000 = 1.8531274795532227
p[k,y,x] 111111 = 1.8531274795532227
p[k,y,x] 000000 = 2.864248514175415
p[k,y,x] 111111 = 3.064248561859131
y = 3
x = 6
p0 = 3
p1 = 6
step = [0.12564097 0.1556096 ]
p[k,y,x] 000000 = 3.125640869140625
p[k,y,x] 111111 = 3.25128173828125
p[k,y,x] 000000 = 6.155609607696533
p[k,y,x] 111111 = 6.311219215393066
y = 3
x = 7
p0 = 3
p1 = 7
step = [0.19941339 0.01530699]
p[k,y,x] 000000 = 3.199413299560547
p[k,y,x] 111111 = 3.3988265991210938
p[k,y,x] 000000 = 7.0153069496154785
p[k,y,x] 111111 = 7.030613899230957
y = 3
x = 8
p0 = 3
p1 = 7
step = [0.19941339 0.01530699]
p[k,y,x] 000000 = 3.144273042678833
p[k,y,x] 111111 = 3.34368634223938
p[k,y,x] 000000 = 7.8614888191223145
p[k,y,x] 111111 = 7.876795768737793
y = 3
x = 9
p0 = 3
p1 = 8
step = [ 0.14427295 -0.13851108]
p[k,y,x] 000000 = 3.052521228790283
p[k,y,x] 111111 = 3.196794271469116
p[k,y,x] 000000 = 8.807019233703613
p[k,y,x] 111111 = 8.668508529663086
y = 4
x = 6
p0 = 3
p1 = 6
step = [0.12564097 0.1556096 ]
p[k,y,x] 000000 = 3.865173101425171
p[k,y,x] 111111 = 3.990813970565796
p[k,y,x] 000000 = 6.147721767425537
p[k,y,x] 111111 = 6.30333137512207
y = 4
x = 7
p0 = 3
p1 = 6
step = [0.12564097 0.1556096 ]
p[k,y,x] 000000 = 3.800049066543579
p[k,y,x] 111111 = 3.925689935684204
p[k,y,x] 000000 = 6.995570182800293
p[k,y,x] 111111 = 7.151179790496826
y = 4
x = 8
p0 = 3
p1 = 7
step = [0.19941339 0.01530699]
p[k,y,x] 000000 = 3.8391151428222656
p[k,y,x] 111111 = 4.0385284423828125
p[k,y,x] 000000 = 7.881188869476318
p[k,y,x] 111111 = 7.896495819091797
y = 4
x = 9
p0 = 3
p1 = 8
step = [ 0.14427295 -0.13851108]
p[k,y,x] 000000 = 3.940187692642212
p[k,y,x] 111111 = 4.084460735321045
p[k,y,x] 000000 = 8.80915355682373
p[k,y,x] 111111 = 8.670642852783203
p = [[[0. 0. 0.30913824 0.30144167 0. 0.
0. 0. 0. 0. ]
[1. 1. 1. 1.1391661 1. 1.
1. 1. 1. 1. ]
[2. 2. 1.839529 1.8531275 2. 2.
2. 2. 2. 2. ]
[3. 3. 3. 3. 3. 3.
3.2512817 3.3988266 3.3436863 3.1967943 ]
[4. 4. 4. 4. 4. 4.
3.990814 3.92569 4.0385284 4.0844607 ]
[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.2538376 2.9911673 4. 5.
6. 7. 8. 9. ]
[0. 1. 2. 2.927513 4. 5.
6. 7. 8. 9. ]
[0. 1. 2.3193696 3.0642486 4. 5.
6. 7. 8. 9. ]
[0. 1. 2. 3. 4. 5.
6.311219 7.030614 7.876796 8.668509 ]
[0. 1. 2. 3. 4. 5.
6.3033314 7.15118 7.896496 8.670643 ]
[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. ]]]
可以发现，此时y和x这两个坐标与p0和p1这两个坐标是不同的。p0和p1是y和x经过整型化以后得到的，而网格格点上的热量累加也是针对的这两个整型化以后的格点上的梯度值。
这里面特别需要理解的地方是：p0和p1不是随意的整型化，而是代表了原格点的热量所来源的方向，因为网格格点矩阵p上更新后的值都是累加的带方向的梯度，这个方向的反方向就代表了其热源的来源。以格点(1,3)为例，它经过第一次迭代后变成了(0.9845, 2.8005)，即它累加了垂直方向上由上到下的热流和水平方向上由左到右的热流，它整型化后变成(0, 2)，即指向了它的热流的来源，下一次迭代时就再累加这个格点上的热流梯度。
经过200次迭代以后，最终的网格格点坐标变成了如下：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. ]]]
复原掩膜的工作由如下函数完成：1
mask = get_masks(p, iscell=cp_mask)
下面进入该函数，看它具体干了啥。
构建流场强度直方图
1 | for i in range(dims): |
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
2edges = [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的值为：1
2
3
4
5
6
7
8
9
10[[1. 1. 0. 0. 1. 1. 1. 1. 1. 1.]
[1. 1. 3. 3. 1. 1. 1. 1. 1. 1.]
[1. 1. 0. 0. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 0. 5. 0. 0.]
[1. 1. 1. 1. 1. 1. 0. 3. 0. 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.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
第一个数值代表pflows中的元素其第一个维度（axis=0，这里没说成实际的Y方向，就是从纯数据的角度来看，不跟具体的物理参数弄混）为-0.5到0.5之间（这里就是特指0），且第二个维度为-0.5到0.5之间（就是0）的数的个数，那么看pflows就是(0, 0)这个元素；第二个数值就代表第一个维度为0，且第二个维度为1的数的个数，发现就是(0, 1)这一对；第三个数值就代表第一个维度为0，且第二个维度为2的数的个数，发现没有这样的数，所以第三个数值为0。，另外，比如第二行的3代表当第一个维度为0.5到1.5之间（就是1）时，且第二个维度为1.5到2.5之间（就是2）的数有三个，分别在pflows的第3个、第13个、第23个元素对。
以上是关于那个直方图矩阵的纯数学理解，对于其实际物理意义，可以想象成这样：如果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.]]]
那么上面的统计矩阵h上的元素应该都为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. 0. 0. 1. 1. 1. 1. 1. 1.]
[1. 1. 3. 3. 1. 1. 1. 1. 1. 1.]
[1. 1. 0. 0. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 0. 5. 0. 0.]
[1. 1. 1. 1. 1. 1. 0. 3. 0. 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.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
看pflows，先观察其第一个维度，对于原来都是0的标识，此时有两个0变成了两个1，而同时这两个标识对应的第二个维度没有变化，表示这两个位置上的米粒由原来0号标识移到了1号标识上，因为pflows的第一个维度实际是垂直方向，所以，就是将这两个位置上的米粒向下移动了一格，显示在直方图统计上，就是这两个位置统计数目为0，而其下的对应位置分别加1；实际上，在同样的垂直位置上，有两个原来标识为2的米粒也变成了标识为1，即这两个米粒往上移动一格，对应于直方图统计就是这两个位置统计为0，而其上对应位置分别加1，所以最终直方图统计上第二行第三列和第四列的统计数目都为3。
对于米粒移动的理解，可以参考如下示意图：
上面是以相对位置角度来理解，另一个角度来看，从绝对位置来看，pflows中的要移动米粒的位置上面的值正是热源所在坐标位置：即要移动到的位置，就是要寻找的热源位置，即相对于原网格点阵有变化的格点上的新坐标就是标明了热源位置，比如，那两个由0变为1的格点，其原坐标分别是(0, 2)和(0, 3)，现在分别变成为(1, 2)和(1, 3)，即它俩指向的热源就在(1, 2)和(1, 3)上。这一点非常重要，后面会用到。
关于直方图统计这一部分理解结束！
寻找局部极大值种子
得到直方图统计结果后，对其进行一个核为5的最大化池化（池化的目的是为了滤波，去掉噪声）：1
2for 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))
接下来是两层循环。 以下代码中的注释和输出是以单次外层循环下的单次内层循环为例。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# 外层循环，目的是为了重复进行多次内层循环，5也是超参数
for iter in range(5):
# 对pix的每个元素进行循环，再次提醒，pix是局部极大值的坐标索引
for k in range(len(pix)):
# 如果是第一次外循环，将每个pix元素由np.array改成list
if iter==0:
pix[k] = list(pix[k])
newpix = []
# 对expand的两个维度进行循环
for i,e in enumerate(expand):
# 第一项是对expand的每个维度中的元素加上一个维度
# 第二项是对pix的当前元素的每个维度中的元素加上一个维度
# 第三项是减去1
epix = e[:,np.newaxis] + np.expand_dims(pix[k][i], 0) - 1
epix = epix.flatten()
# newpix是在两个维度上对epix进行组合
newpix.append(epix)
# 总体来说，newpix就是以pix中的元素为中心的八邻域中的九个邻居像素（加上了它自己）的坐标索引
# 比如pix中某个像素是[1, 2]，那么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)]
newpix = tuple(newpix)
# 然后看之前的直方图统计矩阵h在这九个像素点上是否大于2
# 如果大于2，即代表有热源流入，将igood设为true
# 以pix中的第一个元素[1, 2]为例，其igood为：
# igood = [False False False False True True False False False]
# 这代表了[1, 2]和[1, 3]位置上的热流强度大于2
igood = h[newpix]>2
# 将有热流流入的这个flag带入newpix中，然后再重新存入pix
for i in range(dims):
pix[k][i] = newpix[i][igood]
# 接上例，此时pix变成：
# [[array([1, 1], dtype=int64), array([2, 3], dtype=int64)], array([1, 3], dtype=int64), array([3, 7], dtype=int64)]
# 对比pix之前的数值：
# [array([1, 2], dtype=int64), array([1, 3], dtype=int64), array([3, 7], dtype=int64)]
# 发现就是将pix的第一个元素由原来的局部极大值索引[1, 2]，变成了以它为中心的九个像素中直方图统计大于2的像素坐标索引，注意[[array([1, 1], dtype=int64), array([2, 3], dtype=int64)]]的意思是两个x坐标放在一起，两个y坐标放在一起，实际索引要竖过来看，即[1, 2]和[1, 3]
if iter==4:
pix[k] = tuple(pix[k])
对比一下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]这四个元素（注意竖过来看坐标组合）。
咋一看，发现pix与没有漫水填充前的seeds种子点差不多，即只差了[4, 7]这一个位置。但实际上差别还挺大：
如果只用之前的seeds作为下面的标记点，那么会导致[4, 7]这个位置所接收的热流的来源不会被标记，反映到结果上就会表现出这样：
- 很多像素可能被忽略，没有被正确标记；
- 原本一整块物体可能会被分割得支离破碎，毕竟没法精确控制物体内的热流流向。一个例子，比如一个物体内，发现中心点有两个最终的热流汇聚点，但极大值种子点seeds只有一个，那么如果没有上面的漫水填充后，则只有一个中心点所接收的热流流入被统计；另一个势均力敌的汇聚点所辐射的范围却没有被计入。
标记辐射范围
然后对这些元素进行label：1
2
3M = 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]]
这里也非常巧妙，可以发现：
- 属于同一极值种子点辐射范围的像素，其标记是一致的；
- 不同极值种子点如果其辐射范围是一样的，也不会改变结果，只是会让标记值加1，不影响结果，后面会通过
np.unique(M0, return_inverse=True)进行处理； - 这里就直观看出了漫水填充的作用，它会不只着眼于极大值种子点，还能找到该极大值种子点辐射范围内的其他热流流入点，保证所有向它们流入的像素都能被完全统计。
从上面已经知道，pflows中异常值就代表了热源位置，所以这里根据pflows的索引，得到其上的label值（即将pflows传入M）：1
2
3
4for 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
7class 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
22class 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
7if 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的整合，即对这两个模型的统一调用。
