手撕三维点云中的降采样最远距离采样（farthest point sample）FPS代码

最近在阅读SC3K: Self-supervised and Coherent 3D Keypoints Estimation from Rotated, Noisy, and Decimated Point Cloud Data这篇在三维点云内基于无监督关键点检测的论文并且在阅读源码，顺带记录下对于最远距离采样方法的实现。

一、最远距离采样方法的基本思想

最远距离采样的思想比较简单，直接展开：首先我们会有N个输入的点云，目的是要在这坨点云内通过最远距离采样直接找着n个最远距离点，这个找的过程就是第一次随机确定N个点内的一个点视为第一个最远采样点a，然后其余N-1个点分别与第一个最远采样点计算他们之间的距离，其中在N-1个点中与第一个最远采样点的距离最大的点b就被视为第二个最远距离点，那么此时我们就采集到了两个最远采样点。这时我们整理下，我们要从剩下的N-2个点内（此时还需要找n-2个最远采样点）再次循环找到第三个最远采样点直到找着第n个最远采样点，这里需要注意：我们这里有两个最远采样点了，剩下N-2个点都要与这两个最远采样点a、b分别计算他们之间的距离，这里我们的算法是取距离a、b点间两个距离间的最小值，用这个最小值来代替N-2个点各自与最远点的距离度量，按照这个规则就能在取到第三个最远采样点直到n个最远采样点了~

二、代码实现

def farthest_point_sample(point, npoint): #point表示我们要研究的点云对象，npoint表示要采样的最远采样点个数
    N, D = point.shape
    xyz = point[:,:3] 
    centroids = np.zeros((npoint,))
    distance = np.ones((N,)) * 1e10 
    farthest = np.random.randint(0, N) #从0-（N-1）
    for i in range(npoint):
        centroids[i] = farthest 
        centroid = xyz[farthest, :] 
        dist = np.sum((xyz - centroid) ** 2, -1) 
        mask = dist < distance 
        distance[mask] = dist[mask] 
        farthest = np.argmax(distance, -1) 
    point = point[centroids.astype(np.int32)]
    return point

现在来逐行剖析我在第一次读该代码遇到的疑惑。

N, D = point.shape

这里是解包操作，通过.shape返回研究点云对应的形状，这里输入的点云数据集是N*D维的，N为点云的个数，D为点云的特征数

xyz = point[:, :3]

借助xyz单独提取出输入点云的所有数据，并且用xyz仅提取出输入点云的空间坐标特征，即x轴、y轴、z轴特征，这是专门用来后续计算最远点采样时计算距离的。

centroids = np.zeros((npoint,))

创建n个全0的numpy数组，用来存储后续最远采样点在输入点云的坐标索引

distance = np.ones((N, )) * 1e10

这是用来生成N个全为1*1e10数值的numpy数组，他是用来存储N个点在循环中与当前采样到的所有最远采样点间最小距离的数值，注意是每个点与多个最远采样点间计算的距离，取这个计算的多个距离的最小值存储到distance内，用这个最小距离来通过最大值去提取出最远采样点。

farthest = np.random.randint(0,N) #从0-N-1

这是随机生成0~N-1中一个整数的命令，就将其视为初始化随机生成最远距离点的下标索引即可，从上述整体的程序即可得：该代码只使用了一次，后面在for循环内farthest又被直接重新定义了，这叫python的变量重定义机制，是合理的。

现在进入变量定义好了，那么该进入从0-n-1个最远采样点的循环内噜。

centroids[i] = farthest

还记得前面对centroids的定义吗？他是n个数值为0的numpy数组，farthest是从0-（N-1）的随机生成的一个整数，在i=0时代表生成第一个最远采样点，其中farthest为将其从N个点云的某个点的索引赋予给centroids[0]。

centroid = xyz[farthest, :]

这里定义了一个新的临时变量名centroid，它是用来从输入的点云内提取出最新的一个最远采样点的坐标信息，注意，for循环内每次循环都会计算所有点与最新的一个最远点坐标计算一次距离，并且更新到distance，意味着后面的distance会在循环中不断的更新，并且其对应索引值的数值-距离，对最新最远点间距离会与前面所有先前计算的最远点的距离进行比较，取最小值来不断替换，我一开始一直以为该程序每次循环都会整体性的更新distance，后面发现他每次循环只会以一个最新的最远点坐标来计算距离，以此来逐步更新distance来极大降低计算量！

dist = np.sum((xyz - centroid) ** 2, -1)

注意这里centroid是xyz[farthest, :],所以他的维度是1*3，而xyz的维度为N*3，所以这里用到了python的广播机制，会将centroid扩展到N*3，N是都相同的数值，是计算的欧氏距离，-1代表按行跨列操作，所以注意dist是一个N*1的numpy数组。这里按行跨列操作不太清楚的可以看看我这篇小笔记（python按行跨列解释的小笔记）。通过该代码就能计算出每次采样到新的最远采样点对应的N个距离，这里大家得知道为什么是N个距离：因为xyz也包括了最远采样点！！！所以每次循环的话，我们从N个点内在挑选其他暂时不是最远采样点时也计算了最远采样点对其自身以及对应其他最远采样点的距离！！！只是由于我们在计算每个点与各个最远采样点间的距离最小值来表征该点是否能作为最远采样点的直接信息，由于这个原因就能得知前辈为啥以距离最小值为度量能否判别为最远采样点的信息而不是拿最大值了。

mask = dist < distance

这里最开始的distance是全为1e10的N维数组，由于distance很大，所以dist对应的数值都满足条件，所以mask是全1的N*1数值（注意mask这里是布尔类型数组，放入到后面的代码会返回为0-1值）。那么在第一次循环完成后，这个mask就在最新最远采样点中通过dist<distance这个条件来不断更新distance来作为判断某些点是否为新的最远采样点了！并且这个代码算是整个循环内的点睛之笔了！！！通过这个循环比较，我们就能在循环中每次只与一个新的关键点计算距离来判别最远采样点，原因是通过disk < distance,其中的distance是隐式地将每轮的最远点信息保留住，并且dist会针对最新的最远点来计算他们之间的距离，进而将当前的最远点与以前的最远点进行隐式地比较，来减少大量复杂代码的执行。

distance[mask] = dist[mask]

在第一次循环中distance内的值会全部由1e10变为N个点与第一个随机选取的最远距离点的各自距离。而在第二次循环以及更后面的循环中，新的最远采样点不断生成。这里以某个阶段为例说明：在当前第x个最新最远采样点钟，这个mask会在当此的循环中与之前x-1个最远采样点的最小距离通过dist < distance比对，如果当前第x个最新的采样点与未被划分为最远采样点的其他点计算出的距离 比 先前未被划分为最远采样点计算出的距离小，那么借助mask就会更新，如果大的话mask就会为0，那么distance在对应的位置索引就不会被更新！

farthest = np.argmax(distance, -1)

这里的distance是一维数组，所以-1和0轴索引是一致的，np.argmax为返回distance的最大值的索引值，这样就提取出了最新的最远采样点。

通过这一整套机制，一直循环到n个最远采样点，直至结束。要是哪里没说明白大家可以讨论下哈。

跳出循环后的代码就简单了

point = point[centroids.astype(np.int32)]

centroids表示在N个点云数据中所有最远采样点的位置索引，astype是将该位置索引转变为np.int32数据类型，并最终返回point内所有最远采样点的位置坐标。

总结：一个简单的代码，运用到了python内numpy数组等各种数据结构、轴的使用、在for循环内的变量重定义以及全局变量与临时变量的生命周期、算法简洁性（避免空间复杂度）的设计等，得好好吸收，今天把工位门口的小橘猫喂饱了~。