python - 如何从 Python 中的点列表创建一个以椭圆体为中心的二进制 3 维矩阵？

Question

目标

我有一个长度为 M 的 3 维坐标点数组作为 float 。我想创建一个预定义形状的 3 维 numpy 数组，其中填充了以这些点为中心的给定浮点半径的椭圆体。因为这是用于图像处理的，所以我将数组中的每个值称为“像素”。如果这些椭圆体重叠，我想通过欧氏距离将像素分配到更近的中心。最终输出将是一个 numpy 数组，背景为 0，椭圆体内的像素编号为 1、2、... M，对应于初始坐标列表，类似于 scipy 的 ndimage.label(...) 的输出.

下面，我有一个天真的方法，它考虑输出数组中的每个位置并将其与每个定义的中心进行比较，为任何椭球内的任何像素创建一个值为 1 的二进制数组。然后它使用 scikit-image 来分水岭这个二进制数组。虽然这段代码有效，但对我来说速度太慢了，这既是因为它考虑了每个像素和中心对，也是因为它单独执行分水岭。我怎样才能加快这段代码的速度？

简单的例子

def define_centromeres(template_image, centers_of_mass, xradius = 4.5, yradius = 4.5, zradius = 3.5):
    """ Creates a binary N-dimensional numpy array of ellipsoids.

    :param template_image: An N-dimensional numpy array of the same shape as the output array.
    :param centers_of_mass: A list of lists of N floats defining the centers of spots.
    :param zradius: A float defining the radius in pixels in the z direction of the ellipsoids.
    :param xradius: A float defining the radius in pixels in the x direction of the ellipsoids.
    :param yradius: A float defining the radius in pixels in the y direction of the ellipsoids.
    :return: A binary N-dimensional numpy array.
    """
    out = np.full_like(template_image, 0, dtype=int)
    for idx, val in np.ndenumerate(template_image):
        z, x, y = idx
        for point in centers_of_mass:
            pz, px, py = point[0], point[1], point[2]
            if (((z - pz)/zradius)**2 + ((x - px)/xradius)**2 + ((y - py)/yradius)**2) <= 1:
                out[z, x, y] = 1
                break
    return out

Scikit-image的分水岭函数；通过更改此方法中的代码不太可能找到加速:

def watershed_image(binary_input_image):
    bg_distance = ndi.distance_transform_edt(binary_input_image,
                                             return_distances=True,
                                             return_indices=False)
    local_maxima = peak_local_max(bg_distance, min_distance=1, labels=binary_input_image)
    bg_mask = np.zeros(bg_distance.shape, dtype=bool)
    bg_mask[tuple(local_maxima.T)] = True
    marks, _ = ndi.label(bg_mask)
    output_watershed = watershed(-bg_distance, marks, mask=binary_input_image)
    return output_watershed

小规模示例数据:

zdim, xdim, ydim = 15, 100, 100
example_shape = np.zeros((zdim,xdim,ydim))
example_points = np.random.random_sample(size=(10,3))*np.array([zdim,xdim,ydim])
center_spots_image = define_centromeres(example_shape, example_points)
watershed_spots = watershed_image(center_spots_image)

输出:

center_spots_image, max projected to 2d

watershed_spots, max projected to 2d

请注意，这些图像只是最终 3d 输出数组的 2d 表示。

附加说明

输出数组的典型大小为 31x512x512，即总共 8.1e6 个值，输入坐标的典型大小为 40 个 3 维坐标点。我想针对此秤优化此过程的速度。

我在这个项目中使用了 numpy、scipy 和 scikit-image，我必须坚持使用这些和其他维护良好且记录良好的包。

对于上述代码的可访问性错误或我的解释不够清晰，我们深表歉意。我是一名研究科学家，几乎没有接受过正规的计算机科学培训。

Answer 1

无论是算法还是它的实现确实都留下了一些改进的空间。 @Eric Johnson 已经介绍了后者，现在请允许我通过使用更好的算法来演示该示例的另一个 20 倍加速。

改进:1.在屏蔽之前限制为一个易于计算的边界框。2. 对于重叠分辨率，回收已经为椭球计算完成的距离计算。

代码(假设 Eric 的函数已经定义):

import numpy as np

def rasterise(template_image, centers_of_mass, xradius=4.5, yradius=4.5, zradius=3.5):
    """Creates a labeled N-dimensional numpy array of ellipsoids.

    :param template_image: An N-dimensional numpy array of the same shape as the output array.
    :param centers_of_mass: A list of lists of N floats defining the centers of spots.
    :param zradius: A float defining the radius in pixels in the z direction of the ellipsoids.
    :param xradius: A float defining the radius in pixels in the x direction of the ellipsoids.
    :param yradius: A float defining the radius in pixels in the y direction of the ellipsoids.
    :return: An N-dimensional numpy array, with label `n` for the ellipsoid at index `n-1`.
    """
    sh = template_image.shape
    out = np.zeros(sh,int)
    aux = np.zeros(sh)
    radii = np.array([zradius,xradius,yradius])
    for j,com in enumerate(centers_of_mass,1):
        bboxl = np.floor(com-radii).clip(0,None).astype(int)
        bboxh = (np.ceil(com+radii)+1).clip(None,sh).astype(int)
        roi = out[tuple(map(slice,bboxl,bboxh))]
        roiaux = aux[tuple(map(slice,bboxl,bboxh))]
        logrid = *map(np.square,np.ogrid[tuple(
            map(slice,(bboxl-com)/radii,(bboxh-com-1)/radii,1j*(bboxh-bboxl)))]),
        dst = (1-sum(logrid)).clip(0,None)
        mask = dst>roiaux
        roi[mask] = j
        np.copyto(roiaux,dst,where=mask)
    return out
        

    
    
zdim, xdim, ydim = 15, 100, 100
example_shape = np.zeros((zdim,xdim,ydim))
example_points = np.random.random_sample(size=(10,3))*np.array([zdim,xdim,ydim])
center_spots_image = define_centromeres_labeled(example_shape, example_points)
csi = rasterise(example_shape, example_points)
print("number of pixels dfferent",np.count_nonzero(csi != center_spots_image),"out of",csi.size)
from timeit import timeit
print("Eric",timeit(lambda:define_centromeres_labeled(example_shape, example_points),number=10))
print("loopy",timeit(lambda:rasterise(example_shape, example_points),number=10))

样本运行:

number of pixels dfferent 0 out of 150000
Eric 0.37984768400201574
loopy 0.019632569048553705

警告:

Eric 的代码和我的代码的重叠分辨率略有不同。例如:

不同之处在于(我认为)埃里克(上图)使用标准欧几里德度量，而我(下图)使用椭圆体建议的度量，主要是出于投机取巧，但也因为这甚至可能是正确的做法。将其切换是可能的，但会降低速度。