需要一种从大量 3D 坐标中绘制平面的有效方法

Posted

技术标签:

【中文标题】需要一种从大量 3D 坐标中绘制平面的有效方法【英文标题】:Need an efficient way to plot planes from large sets of 3D coordinates 【发布时间】:2021-07-25 01:08:10 【问题描述】:

我在 2D 相机上收集了一些探测器数据,然后将其转换为实验室框架,因此我最终得到图像中每个像素的 (x^2+y^2) 和 z 坐标。但是随后对象围绕它正常旋转,并且每次旋转都有一个 img 。我将旋转矩阵应用于 (x^2+y^2) 以获得每个 imgxy 矩阵,所以我最终为每个图像/角度得到类似的结果。所以每个像素都有一个 3D 位置和强度。

z                  x            y          img
444444444     123456789     123456789    123456789                  
333333333     123456789     123456789    423466789
222222222     123456789     123456789    223256789
111111111     123456789     123456789    523456689

然后我想做的是提取一个平面,即为给定的 z 范围绘制 x、y 的地图。

以下问题稍微复杂一些:

labframe 实际上是弯曲的,所以我不能依赖 x 和 y 的每一行都相同。 图片大小约为 2048x2048x32bits (Tiff) - 可以有 1000 张图片。

我当前的解决方案是使用 CUDA/Numba,我有一个函数可以计算给定角度的z,x,y,img,所以我对所有角度都这样做。每次我然后切片一些行,并用xyimg 值扩展一个列表。然后使用scipy.interpolate.griddata 给出一个二维地图。 griddata 也很慢,GPU 上的任何东西都可能会更好。

整个过程很慢,所以我正在寻找更好的解决方案,或者图书馆可能已经这样做了? CUDA 代码看起来像这样,它本身并没有那么慢:

#constants are q0, angi, rot_direction, SDD, k0, Binv
@cuda.jit
    def detector_to_hkl_kernel(h_glob,k_glob,l_glob,omega_rad):
        #get the current thread position
        j,i = cuda.grid(2)

        if j < h_glob.shape[0] and i < h_glob.shape[1]:
            delta_z= (q0[1]-j)*pixel_y  #real-space dinstance from centre pixel y
            delta_x = (i-q0[0])*pixel_x  #real-space dinstance from centre pixel x
            delR = math.sqrt(delta_x**2 + delta_z**2)            
            dist = math.sqrt(delta_x**2+SDD**2 + delta_z**2) #distance to pixel      

            #lab coorindates of pixel in azimuthal angles
            del_pix  = math.atan(delta_x/ SDD)
            gam_pix = math.atan(delta_z/math.sqrt(delta_x**2 + SDD**2))-angi*math.cos(del_pix)
                            
            #lab coordinates in momenturm transfer                                  
            qx = k0*(math.cos(gam_pix)*math.cos(del_pix)-math.cos(angi))
            qy = k0*(math.cos(gam_pix)*math.sin(del_pix)) 
            qz = k0*(math.sin(gam_pix)+math.sin(angi))

            so = math.sin(rotDirection*omega_rad)
            co = math.cos(rotDirection*omega_rad)
            # we deal with the angle of incidence in the momentum transfer calc
            # so that part of the rotation matrix can be fixed
            ci = 1 #math.cos(angi) 
            si = 0 #math.sin(angi)

            #rotation matrix
            hphi_1 = so*(ci*qy+si*qz)+co*qx
            hphi_2 = co*(ci*qy+si*qz)-so*qx
            hphi_3 = ci*qz-si*qy
                
            #H= Binv dot Hphi 
            # compute the dot product manually 
            h_glob[j,i] = Binv[0][0]*hphi_1+Binv[0][1]*hphi_2+Binv[0][2]*hphi_3
            k_glob[j,i] = Binv[1][0]*hphi_1+Binv[1][1]*hphi_2+Binv[1][2]*hphi_3
            l_glob[j,i] = Binv[2][0]*hphi_1+Binv[2][1]*hphi_2+Binv[2][2]*hphi_3              
            
        
    h_global_mem  = cuda.to_device(np.zeros((pixel_count_y,pixel_count_x)))
    k_global_mem  = cuda.to_device(np.zeros((pixel_count_y,pixel_count_x)))
    l_global_mem  = cuda.to_device(np.zeros((pixel_count_y,pixel_count_x)))                  

    # Configure the blocks
    threadsperblock = (16, 16)
    blockspergrid_x = int(math.ceil(pixel_count_y / threadsperblock[0]))
    blockspergrid_y = int(math.ceil(pixel_count_x / threadsperblock[1]))
    blockspergrid = (blockspergrid_x, blockspergrid_y)

    detector_to_hkl_kernel[blockspergrid, threadsperblock](h_global_mem,k_global_mem,l_global_mem, omega_rad)        
    return [h_global_mem.copy_to_host(),k_global_mem.copy_to_host(),l_global_mem.copy_to_host()]  

【问题讨论】:

代码对我来说似乎并不算太​​糟糕(除了可能不可避免的三角函数调用)。您使用什么 GPU? 我看不到内核中在哪里使用了delRdist。我们希望内核将它们优化掉,但是用冗余代码来诱惑命运有什么意义呢? 啊,是的,我错过了当我为发布进行简化时,还有另一个强度校正数组,它是针对它们用于的每个像素计算的。是的,主要问题实际上是对大量图像执行此操作,然后从最终结果中取出飞机。 使用 GTX 1660 super,但只是因为我的 3080 上周爆炸了(我认为这无关哈哈)。我的 quadro 现在都 4 或 5 岁了。 【参考方案1】:

首先,请注意您在此处使用双精度,而主流的中端消费类 GPU 非常慢 来计算双精度浮点数。事实上,GTX 1660 Super GPU 的计算能力为 5027 GFlops 的简单精度和只有 157 GFlops 的双精度(慢 32 倍)。一种简单的解决方案是在代码中使用简单精度浮点数,方法是指定dtype=np.float32 或使用array.astype(np.float32) 转换数组。如果您不能使用简单精度或混合精度,另一种昂贵的解决方案可能是使用专用的专业 GPU。

此外,可以提前预先计算几个表达式并存储在常量中。这包括例如math.cos(angi)math.sin(angi)1.0/SDD。其他一些表达式可以存储在临时变量中,因为编译器可能无法有效地分解代码(主要是因为trigonometric functions)。

此外,三角函数通常非常昂贵,尤其是当您希望计算符合 IEEE-754 时(math.xxx 调用可能就是这种情况)。您可以改用近似值。 CUDA 提供了 __cosf__sinf__tanf 内在函数,它们应该更快(但如果使用它们,请注意结果)。我不确定你是否可以直接调用它们,但你可以将参数 fastmath=True 添加到 JIT 装饰器中,它可以为你做到这一点。

我认为使用 32x8 的 2D 线程块可能会更快一些,因为线程封装在包含 32 个线程和 GPU 的 warp 中。但最好的解决方案是检查许多不同块大小的性能。

如果所有这些还不够,您可以尝试使用共享内存来减少每个 bloc 执行的指令量,因为每个 bloc 会多次重新计算某些表达式。

【讨论】:

我注意到你再次发生了 OP 在得到答案后删除了问题 ***.com/questions/67359104/… 。我认为在这种情况下,您可以标记问题。尽管如此,我还是投了反对票 @dreamcrash 谢谢!我没有注意到这一点。我标记了这个问题,因为这确实有点可疑。很快,我就能看到已删除的帖子了:)。 不知道你是什么意思删除问题.. 也感谢 Jerome 的提示,我一定会实现它们。我有点忙于其他工作,所以还没有机会回来做这个。但主要瓶颈仍在从此类 3D 数据中提取平面/切片。我想我也需要编写一些代码来在 GPU 上进行直方图/插值 - 如果我能并行化自己的话。 最后我通过合并到网格并使用原子添加解决了这个问题

以上是关于需要一种从大量 3D 坐标中绘制平面的有效方法的主要内容,如果未能解决你的问题,请参考以下文章

将 3D 光标捕捉到平面的不透明部分(搅拌机)

给定表面法线,找到 3D 平面的旋转

查找/重新映射 OpenGL ES 坐标平面的边界

leap motion怎么对应屏幕上的位置

如何使用4个角的位置获得3d平面的角度

通过 3D 向量拟合线并找到与平面的交点