python - 循环遍历大型 .tif 堆栈(图像光栅)并提取位置

Question

我想运行一个大的tif堆栈+1500帧并提取每帧的局部最大值的坐标。下面的代码可以完成这项工作，但是对于大文件来说速度非常慢。当在较小的位(例如 20 帧)上运行时，每个帧几乎立即完成 - 当在整个数据集上运行时，每个帧需要几秒钟。

有什么解决方案可以运行更快的代码吗？我认为这是由于加载了大的 tiff 文件 - 但最初应该只需要一次？

我有以下代码:

from pims import ImageSequence
from skimage.feature import peak_local_max

def cmask(index,array):
  radius = 3
  a,b = index
  nx,ny = array.shape
  y,x = np.ogrid[-a:nx-a,-b:ny-b]
  mask = x*x + y*y <= radius*radius

  return(sum(array[mask])) # number of pixels

images = ImageSequence('tryhard_red_small.tif')


frame_list = []
x = []
y = []
int_liposome = []
BG_liposome = []

for i in range(len(images[0])):
    tmp_frame = images[0][i]

    xy = pd.DataFrame(peak_local_max(tmp_frame, min_distance=8,threshold_abs=3000))
    x.extend(xy[0].tolist())
    y.extend(xy[1].tolist())

    for j in range(len(xy)):
        index = x[j],y[j]    
        int_liposome.append(cmask(index,tmp_frame))

    frame_list.extend([i]*len(xy))
    print "Frame: ", i, "of ",len(images[0])

features = pd.DataFrame(
        {'lip_int':int_liposome,
         'y'  : y,
         'x'  : x,
         'frame'  : frame_list})

Answer 1

您是否尝试过对代码进行分析，例如在 ipython 中使用 %prun 或 %lprun ？这会准确地告诉您速度下降的位置。

如果没有 tif 堆栈，我无法制作自己的版本，但我怀疑问题在于您正在使用列表来存储所有内容。每次执行追加或扩展时，Python 都必须分配更多内存。您可以尝试首先获取最大值的总数，然后分配输出数组，然后重新运行以填充数组。类似下面的内容

# run through once to get the count of local maxima
npeaks = (len(peak_local_max(f, min_distance=8, threshold_abs=3000))
          for f in images[0])
total_peaks = sum(npeaks)

# allocate storage arrays and rerun
x = np.zeros(total_peaks, np.float)
y = np.zeros_like(x)
int_liposome = np.zeros_like(x)
BG_liposome = np.zeros_like(x)

frame_list = np.zeros(total_peaks, np.int)
index_0 = 0
for frame_ind, tmp_frame in enumerate(images[0]):
    peaks = pd.DataFrame(peak_local_max(tmp_frame, min_distance=8,threshold_abs=3000))
    index_1 = index_0 + len(peaks)
    # copy the data from the DataFrame's underlying numpy array
    x[index_0:index_1] = peaks[0].values
    y[index_0:index_1] = peaks[1].values
    for i, peak in enumerate(peaks, index_0):
        int_liposome[i] = cmask(peak, tmp_frame)
    frame_list[index_0:index_1] = frame_ind
    # update the starting index
    index_0 = index_1
    print "Frame: ", frame_ind, "of ",len(images[0])