加快应用于 numpy 3-d 的 scipy ndimage 测量

Question

我有多个大型标记的 numpy 二维数组 (10 000x10 000)。对于每个标签（具有相同编号的连接单元格），我想根据另一个 numpy 3-d 数组（平均值、标准、最大值等）的值计算多个测量值。如果将 3-d numpy 转换为 2-d，则使用 scipy.ndimage.labeled_comprehension 工具可以做到这一点。然而，由于标签的数量和数组的大小相当大，计算需要相当长的时间。我当前的代码似乎是多余的，因为我现在正在为输入图像的每个 3 维迭代相同的标签。我想知道是否有办法加快我的代码速度（例如，将三个 scipy.ndimage.labeled_comprehension 计算合并为一个计算）。

使用形状为 (4200,3000,3) 和 283047 个标签的测试数据集，计算耗时 10:34 分钟

---

测试数据

example_labels=np.array([[1, 1, 3, 3],
   [1, 2, 2, 3],
   [2, 2, 4, 4],
   [5, 5, 5, 4]])

unique_labels=np.unique(example_labels)
value_array=np.arange(48).reshape(4,4,3)

当前代码和所需输出

def mean_std_measurement(x):
    xmean = x.mean()
    xstd = x.std()
    vals.append([xmean,xstd])

def calculate_measurements(labels, unique_labels, value_arr):
    global vals
    vals=[]
    ndimage.labeled_comprehension(value_array[:,:,0],labels,unique_labels,mean_std_measurement,float,-1)
    val1=np.array(vals)
    vals=[]
    ndimage.labeled_comprehension(value_array[:,:,1],labels,unique_labels,mean_std_measurement,float,-1)
    val2=np.array(vals)
    vals=[]
    ndimage.labeled_comprehension(value_array[:,:,2],labels,unique_labels,mean_std_measurement,float,-1)
    val3=np.array(vals)
    return np.column_stack((unique_labels,val1,val2,val3))

>>> print calculate_measurements(example_labels,unique_labels,value_array)
array([[  1.        ,   5.        ,   5.09901951,   6.        ,
      5.09901951,   7.        ,   5.09901951],
   [  2.        ,  21.        ,   4.74341649,  22.        ,
      4.74341649,  23.        ,   4.74341649],
   [  3.        ,  12.        ,   6.4807407 ,  13.        ,
      6.4807407 ,  14.        ,   6.4807407 ],
   [  4.        ,  36.        ,   6.4807407 ,  37.        ,
      6.4807407 ,  38.        ,   6.4807407 ],
   [  5.        ,  39.        ,   2.44948974,  40.        ,
      2.44948974,  41.        ,   2.44948974]])

Answer 1

mean_std_measurement 中附加的列表可能会在处理大型数组时形成瓶颈。 labeled_comprehension 将自行返回一个数组，因此第一步就是让 scipy 在后台处理数组构造。

labeled_comprehension 只能应用输出单个值的函数——我怀疑这就是你首先使用列表构造的原因——但我们可以通过让函数输出一个复杂值来欺骗这一点。另一种选择是使用结构化数据类型作为输出，如果我们返回超过 2 个值，这将是必要的。

import numpy as np
from scipy import ndimage
example_labels=np.array([[1, 1, 3, 3],
   [1, 2, 2, 3],
   [2, 2, 4, 4],
   [5, 5, 5, 4]])

unique_labels=np.unique(example_labels)
value_array=np.arange(48).reshape(4,4,3)
# return a complex number to get around labled_comprehension limitations
def mean_std_measurement(x):
    xmean = x.mean()
    xstd = x.std()
    return np.complex(xmean, xstd)

def calculate_measurements(labels, unique_labels, value_array):
    val1 = ndimage.labeled_comprehension(value_array[:,:,0],labels,unique_labels,mean_std_measurement,np.complex,-1)
    val2 = ndimage.labeled_comprehension(value_array[:,:,1],labels,unique_labels,mean_std_measurement,np.complex,-1)
    val3 = ndimage.labeled_comprehension(value_array[:,:,2],labels,unique_labels,mean_std_measurement,np.complex,-1)
    # convert the complex numbers back into reals
    return np.column_stack((np.unique(labels),
                            val1.real,val1.imag,
                            val2.real,val2.imag,
                            val3.real,val3.imag))