【问题标题】:increase speed for looping through numpy array提高循环遍历 numpy 数组的速度
【发布时间】:2019-08-11 06:21:54
【问题描述】:

我正在尝试在地面分类后分割 LiDAR 点云。我正在使用 numpy 创建点云(pc)的“图像”,并循环遍历 numpy 数组。我想加快循环或一起避免它。我将使用图像分割技术,但首先我需要运行这段代码来创建一个“图像”,这是需要一段时间的部分。有没有办法提高这个循环的速度或避免它?


import numpy as np
from math import ceil, floor


'''In this case:
pc = point cloud (X,Y,Z values)'''

# point cloud is in the numpy array, pc
minx,maxx,miny,maxy = floor(np.min(pc[:,0]-1)),ceil(np.max(pc[:,0]+1)),floor(np.min(pc[:,1]-1)),ceil(np.max(pc[:,1]+1))# x,y bounding box

# grid x and y direction (resolution: 0.2 meters)
gridx = np.linspace(minx,maxx,int((maxx - minx+0.2)*5),endpoint=True) 
gridy = np.linspace(miny,maxy,int((maxy - miny +0.2)*5),endpoint=True)

#shape of the new image with 0.2 meter resolution.
imgx,imgy = int((maxx-minx+0.2)*5),int((maxy - miny +0.2)*5)

# this is what will be created at the end.  It will be a binary image.
img = np.zeros((imgx,imgy))

#loop through array to generate image (this is the part that takes a while)
for x,i in enumerate(gridx):
    for y,j in enumerate(gridy):

# Test if there any points in this "grid"
        input_point = pc[np.where(((pc[:,0]>i) & (pc[:,0]<i+1))& ((pc[:,1]>j) & (pc[:,1]<j+1)))]
# if there are points, give pixel value 1.
        if input_point.shape[0]!=0:
            img[x,y]=1

print('Image made')


谢谢。

【问题讨论】:

    标签: python arrays loops numpy lidar


    【解决方案1】:

    这是一个矢量化版本,它在随机测试集上产生相同的输出:

    import numpy as np
    from math import ceil, floor
    import time
    
    width = 0.2
    
    t = [time.time()]
    
    pc = np.random.uniform(-10, 10, (100, 3))
    
    # point cloud is in the numpy array, pc
    minx,maxx,miny,maxy = floor(np.min(pc[:,0]-1)),ceil(np.max(pc[:,0]+1)),floor(np.min(pc[:,1]-1)),ceil(np.max(pc[:,1]+1))# x,y bounding box
    
    # grid x and y direction (resolution: 0.2 meters)
    gridx = np.linspace(minx,maxx,int((maxx - minx+0.2)*5),endpoint=True) 
    gridy = np.linspace(miny,maxy,int((maxy - miny +0.2)*5),endpoint=True)
    
    #shape of the new image with 0.2 meter resolution.
    imgx,imgy = int((maxx-minx+0.2)*5),int((maxy - miny +0.2)*5)
    
    print('Shared ops done')
    t.append(time.time())
    
    # this is what will be created at the end.  It will be a binary image.
    img = np.zeros((imgx,imgy))
    
    #loop through array to generate image (this is the part that takes a while)
    for x,i in enumerate(gridx):
        for y,j in enumerate(gridy):
    
    # Test if there any points in this "grid"
            input_point = pc[np.where(((pc[:,0]>i) & (pc[:,0]<i+width))& ((pc[:,1]>j) & (pc[:,1]<j+width)))]
    # if there are points, give pixel value 1.
            if input_point.shape[0]!=0:
                img[x,y]=1
    
    t.append(time.time())            
    print('Image made')
    
    if width == 0.2:
        img2 = np.zeros((imgx, imgy), 'u1')
        x2, y2 = (((pc[:, :2] - (minx, miny)) * (5, 5))).astype(int).T
        img2[x2, y2] = 1
    
    elif width == 1:
        img2 = np.zeros((imgx+4, imgy+4), 'u1')
        x2, y2 = (((pc[:, :2] - (minx, miny)) * (5, 5))).astype(int).T
    
        np.lib.stride_tricks.as_strided(img2, (imgx, imgy, 5, 5), 2 * img2.strides)[x2, y2] = 1
        img2 = img2[4:, 4:]
    
    t.append(time.time())
    print('Image remade')
    
    print('took', np.diff(t), 'secs respectively')
    assert((img2==img).all())
    print('results equal')
    

    您的代码生成 5x5 像素。这是故意的吗?我必须有点技巧才能重现它。

    更新:添加了一个制作普通像素的版本。

    示例运行:

    Shared ops done
    Image made
    Image remade
    took [2.29120255e-04 1.54510736e-01 1.44481659e-04] secs respectively
    results equal
    

    【讨论】:

    • 谢谢!我试图通过获取最大和最小 y 值并将它们分成 0.2 米来获得 0.2 米的像素。这不是真的吗?抱歉,我试图理解你的意思。
    • @Steven 但是您检查的间隔是 i 到 i+1 和 j 到 j+1,它们在每个方向上跨越 5 个网格单元。
    • 我并没有尝试跨越 5 个网格单元。我试图做的是根据 0.2 间隔进行检查。例如,如果i = 300.0 和j = 400.0,它将检查x 中的300.0 - 300.2 和y 中的400.0 - 400.2 中是否存在点。它不这样做吗?
    • 哦,呵呵!!那不是它的作用。射击......我想我必须改变它。我的意思是间隔中的下一个数字,我想我只是添加 1,这将使它成为 5 个像素,因为每个像素都是 0.2 米。谢谢你接听。
    • @Steven 实际上,这让事情变得更简单了。我已经更新了答案。
    猜你喜欢
    • 1970-01-01
    • 2017-08-21
    • 1970-01-01
    • 2021-12-02
    • 2021-11-07
    • 1970-01-01
    • 2015-05-19
    • 2015-03-24
    • 1970-01-01
    相关资源
    最近更新 更多