无法将回溯算法移植到cuda内核

问题描述

我一直在尝试使用numba.cuda优化我的图像处理库，大部分都是成功的。但是，以下回溯功能使我非常头疼。

该函数接收类型为numpy.ndarray的{{1}}和形状为uint8的图像。它为(r,c,3)类型和形状float32的图像计算能量图。然后回溯整个能量图，计算出函数返回的回溯数组和最小能量数组。

(r,c)

下面是我尝试将此功能移植到cuda内核的尝试。注意import numpy as np from numba import cuda FILTER_DU = np.stack([np.array([ [1.0,2.0,1.0],[0.0,0.0,0.0],[-1.0,-2.0,-1.0],])] * 3,axis=2) FILTER_DV = np.stack([np.array([ [1.0,[2.0,-2.0],[1.0,axis=2) @cuda.jit def b_convolve(result,img): filter_du = cuda.const.array_like(FILTER_DU) filter_dv = cuda.const.array_like(FILTER_DV) # 2D coords of the current thread i,j = cuda.grid(2) # Ignore thread if coords are outside image img_x,img_y,img_z = img.shape if (i >= img_x) or (j >= img_y): return delta_x = filter_du.shape[0] // 2 delta_y = filter_du.shape[1] // 2 #delta_z = filter_du.shape[1] // 2 <- not needed since 3rd dim is always 3 long (R,G,B) # The result at coordinates (i,j) is equal to # abs(sum_{k,l,m} filter_du[k,m] * img[i-k+delta_x,j-l+delta_y,1-delta_z]) + # abs(sum_{k,m} filter_dv[k,1-delta_z]) # with k,l and m going through the whole mask array s1=0 s2=0 for k in range(filter_du.shape[0]): for l in range(filter_du.shape[1]): i_k = i - k + delta_x j_l = j - l + delta_y # Check if (i_k,j_k) coordinates are inside the image: if (0 <= i_k < img_x) and (0 <= j_l < img_y): s1 += filter_du[k,0] * img[i_k,j_l,1] s1 += filter_du[k,1] * img[i_k,0] s1 += filter_du[k,2] * img[i_k,-1] s2 += filter_dv[k,1] s2 += filter_dv[k,0] s2 += filter_dv[k,-1] result[i,j] = abs(s1)+abs(s2) def calc_energy(img): img = np.ascontiguousarray(img.astype('float32')) energy_map = np.empty(img.shape[:2]).astype('float32') blockdim = (32,32) griddim = (img.shape[0] // blockdim[0] + 1,img.shape[1] // blockdim[1] + 1) b_convolve[griddim,blockdim](energy_map,img) return energy_map def create_backtrack(img): r,_ = img.shape energy_map = calc_energy(img) M = energy_map.copy() backtrack = np.zeros_like(M,dtype=np.int) for i in range(1,r): for j in range(0,c): # Handle the left edge of the image,to ensure we don't index a -1 if j == 0: idx = np.argmin(M[i-1,j:j + 2]) backtrack[i,j] = idx + j min_energy = M[i-1,idx + j] else: idx = np.argmin(M[i - 1,j - 1:j + 2]) backtrack[i,j] = idx + j - 1 min_energy = M[i - 1,idx + j - 1] M[i,j] += min_energy return M,backtrack被idx + j取代，因为idx已经占了cuda_argmin（或者至少应该占）。

查看给定4x4 image的函数的输出，我们可以看到结果甚至不接近匹配。

@cuda.jit('int32(float32[:,:],int32,int32)',device=True)
def cuda_argmin(M,i,j,k):
    min_val = M[i,j]
    min_ind = j
    for x in range(j+1,k):
        if M[i,x] < min_val:
            min_val = M[i,x]
            min_ind = x
    return min_ind

@cuda.jit('void(int32[:,float32[:,:])')
def cuda_backtrack(backtrack,M):
    i,j = cuda.grid(2)
    
    if i == 0:
        return
        
    min_energy = 0
    if j == 0:
        idx = cuda_argmin(M,i-1,j+2)
        backtrack[i,j] = idx
        min_energy = M[i-1,idx]
    else:
        idx = cuda_argmin(M,j-1,j] = idx-1
        min_energy = M[i-1,idx-1]
    M[i,j] += min_energy

def create_backtrack(img):
    r,_ = img.shape
    energy_map = calc_energy(img)
    
    M = energy_map.copy()
    backtrack = np.zeros_like(M,dtype=np.int)
    
    blockdim = (32,32)
    griddim = (r // blockdim[0] + 1,c // blockdim[1] + 1)
    cuda_backtrack[griddim,blockdim](backtrack,M)
    
    return M,backtrack

# Regular backtracking function
M = array([[1750.,1622.,1752.,1176.],[3364.,2982.,3124.,2538.],[5060.,4754.,4228.,4212.],[6278.,6528.,6354.,5866.]],dtype=float32)

backtrack = array([[0,0],[1,1,3,3],2,3]])

我显然做错了事，但我不知道在哪里。输出的差异如何？

编辑：根据@talonmies的建议，我在代码中包含了# Cuda backtrack function M = array([[4.2e-45,2.9e-44,3.1e-44,3.2e-44],[3.6e-44,1.5e-44,1.7e-44,2.0e-44],[2.2e-44,2.8e-44,2.9e-44],[3.2e-44,3.8e-44,4.2e-44,4.3e-44]],dtype=float32) backtrack = array([[ 0,-1,[ 4,2],[ 2,4,6],[ 6,9,10,11]])函数，因此只要您拥有calc_energy和numba，它就至少可以运行。已安装。

解决方法

您的代码中至少有3个问题，其中第一个已经指出：

您将M用作输入和输出。 M的输出值取决于M的其他值。由于CUDA不提供有关线程执行顺序的任何语句，因此将CUDA线程分配给M的每个输出值通常不会得到与串行版本相同的结果。但是，我们可以观察到给定行（M）的i值仅取决于前一行（M）的i-1值。这意味着我们可以跨行并行，但不能跨列并行（至少不容易）。
在创建backtrack数组时，您似乎在正确处理左边缘（j==0），但没有正确处理右边缘（j==c-1）。这种编码错误将导致cuda输出和numpy输出之间的数值差异。
您的两个argmin函数不会返回相同的内容，因此您没有正确地说明这一点。 numpy argmin返回0、1或2。向其中添加j-1时，会得到j-1，j或j+1作为{{1}的可能值}条目。您的cuda argmin函数返回backtrack，j-1或j。减去1时，这些值与numpy不匹配。

如建议的那样，我们可以通过并行执行来解决第一个问题。但是，在继续计算下一行之前，必须对每行进行完全计算。为了有效地使用它，我们需要跨所有CUDA线程的同步点。 CUDA在块级提供了此功能，但是在网格级CUDA中未通过numba CUDA公开此功能。

下面的代码足以解决上述所有问题，因此至少两个不同方法的输出匹配。但是，其余问题与项目1有关。此方法仅使用单个线程块，因此我们可以利用所有线程的每行同步。对于GPU的性能而言，这并不是一个好的最终解决方案。如果您要执行多组这些操作，则可以通过为每个块分配一个图像处理任务来并行执行这些操作，并以此方式获得良好的性能。

j+1

我并不是说上面的代码没有缺陷，也不适合任何特定目的。它主要是您的代码，有一些“已修复”的项目概述了我发现的问题。特别是，如果您想将此代码扩展到上面给出的1-threadblock示例（或宽度大于1024的图像）之外，则有必要进行一些额外的设计工作，

cuda cuda numba python