2

我有以下数据数组,有 200 万个条目:

[20965  1239   296   231    -1    -1 20976  1239   299   314   147   337
   255   348    -1    -1 20978  1239   136   103   241   154    27   293
    -1    -1 20984  1239    39   161   180   184    -1    -1 20990  1239
   291    31   405    50   569   357    -1    -1 20997  1239   502    25
   176   215   360   281    -1    -1 21004  1239    -1    -1 21010  1239
   286   104   248   252    -1    -1 21017  1239   162    38   331   240
   368   363   321   412    -1    -1 21024  1239   428   323    -1    -1
 21030  1239    -1    -1 21037  1239   325    28   353   102   477   189
   366   251   143   452 ... ect

该阵列包含 CCD 芯片上光子的 x,y 坐标,我想通过阵列并将所有这些光子事件添加到尺寸等于 CCD 芯片的矩阵中。

格式如下:number number x0 y0 x1 y1 -1 -1. number我不太关心的两个条目,x0 y0 等。是我想出去的。条目是指示新帧的-1分隔符,在这些之后总是有 2 个“数字”条目。

我制作了这段代码,它确实有效:

i = 2
pixels = np.int32(data_height)*np.int32(data_width)
data = np.zeros(pixels).reshape(data_height, data_width)

while i < len(rdata):
    x = rdata[i]
    y = rdata[i+1]

    if x != -1 and y != -1:
        data[y,x] = data[y,x] + 1
        i = i + 2
    elif x == -1 and y == -1:
        i = i + 4
    else:
        print "something is wrong"
        print i
        print x
        print y

rdata是我的原始数组。data是从零开始的结果矩阵。while 循环从第一个x坐标开始,在索引 2 处,然后如果找到两个连续的-1条目,它将跳过四个条目。

该脚本运行良好,但运行需要 7 秒。我怎样才能加快这个脚本?我是 python 的初学者,从学习 python 最难的方式来看,我知道应该避免 while 循环,但重写为 for 循环更慢!

for i in range(2, len(rdata), 2):

    x = rdata[i]
    y = rdata[i+1]

    if x != -1 and y != -1:
        px = rdata[i-2]
        py = rdata[i-1]

        if px != -1 and py != -1:
            data[y,x] = data[y,x] + 1

也许有人可以想到一种更快的方法,类似于np.argwhere(rdata == -1)并使用此输出来提取xy坐标的位置?


更新:感谢所有答案!

我使用了 askewchan 的方法来保存帧信息,但是,由于我的数据文件有 300000 帧长,当我尝试生成一个尺寸为 (300000、640、480) 的 numpy 数组时出现内存错误。我可以通过制作一个生成器对象来解决这个问题:

def bindata(splits, h, w, data):

    f0=0
    for i,f in enumerate(splits):
        flat_rdata = np.ravel_multi_index(tuple(data[f0:f].T)[::-1], (h, w))
        dataslice = np.zeros((w,h), dtype='h')
        dataslice = np.bincount(flat_rdata, minlength=pixels).reshape(h, w)
        f0 = f
        yield dataslice

然后,我使用 Gohlke 的tifffile.py的修改版本从数组中创建一个 tif,以从数据中生成一个 tiff 文件。它工作正常,但我需要找到一种压缩数据的方法,因为 tiff 文件大于 4gb(此时脚本崩溃)。我有非常稀疏的数组,640*480 全零,每帧几十个,原始数据文件是 4MB,所以应该可以进行一些压缩。

4

3 回答 3

6

听起来您想要做的只是做一些布尔索引魔术来摆脱无效的帧内容,然后当然将像素相加。

rdata = rdata.reshape(-1, 2)
mask = (rdata != -1).all(1)

# remove every x, y pair that is after a pair with a -1.
mask[1:][mask[:-1] == False] = False
# remove first x, y pair
mask[0] = False

rdata = rdata[mask]

# Now need to use bincount, [::-1], since you use data[y,x]:
flat_rdata = np.ravel_multi_index(tuple(rdata.T)[::-1], (data_height, data_width))

res = np.bincount(flat_rdata, minlength=data_height * data_width)
res = res.reshape(data_height, data_width)
于 2013-05-10T15:24:27.503 回答
2

使用它来删除-1s 和numbers:

rdata = np.array("20965  1239   296   231    -1    -1 20976  1239   299   314   147   337 255   348    -1    -1 20978  1239   136   103   241   154    27   293 -1    -1 20984  1239    39   161   180   184    -1    -1 20990  1239 291    31   405    50   569   357    -1    -1 20997  1239   502    25 176   215   360   281    -1    -1 21004  1239    -1    -1 21010  1239 286   104   248   252    -1    -1 21017  1239   162    38   331   240 368   363   321   412    -1    -1 21024  1239   428   323    -1    -1 21030  1239    -1    -1 21037  1239   325    28   353   102   477   189 366   251   143   452".split(), dtype=int)

rdata = rdata.reshape(-1,2)
splits = np.where(np.all(rdata==-1, axis=1))[0]
nonxy = np.hstack((splits,splits+1))
data = np.delete(rdata, nonxy, axis=0)[1:]

现在,使用@seberg 的部分方法将 xy 列表转换为数组,您可以制作一个 3D 数组,其中每个“层”都是一个框架:

nf = splits.size + 1            # number of frames
splits -= 1 + 2*np.arange(nf-1) # account for missing `-1`s and `number`s
datastack = np.zeros((nf,h,w))
f0 = 0                          # f0 = start of the frame
for i,f in enumerate(splits):   # f  = end of the frame
    flat_data = np.ravel_multi_index(tuple(data[f0:f].T)[::-1], (h, w))
    datastack[i] = np.bincount(flat_rdata, minlength=h*w).reshape(h, w)
    f0 = f

现在,是一个显示数据的第 th 帧的datastack[i]2D 数组。i

于 2013-05-10T15:24:33.133 回答
0

如果x0, y0, x1, y1 != -1你能不做类似filter(lambda a: a != -1, rdata)的事情,然后不打扰ifs吗?这可以加快您的代码速度。

于 2013-05-10T15:25:35.463 回答