如何解决需要一种从大量 3D 坐标中绘制平面的有效方法
我在 2D 相机上收集了一些探测器数据，然后将其转换为实验室帧，这样我最终会得到图像中每个像素的 (x^2+y^2) 和 z 坐标。但是然后对象围绕它正常旋转，并且每次旋转都有一个 img 。我将旋转矩阵应用于 (x^2+y^2) 以获得每个 x 的 y 和 img 矩阵，所以我最终会为每个图像/角度得到类似的结果。所以每个像素都有一个 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，所以我这样做这对于所有的角度。每次我然后切片多行，并扩展一个列表，其中包含散乱的 x,img 值。然后使用 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 非常慢计算双精度浮点数。事实上，GTX 1660 Super GPU 的简单精度计算能力为 5027 GFlops，双精度仅为 157 GFlops（慢 32 倍）。一种简单的解决方案是在您的代码中使用简单精度浮点数，方法是指定 if 或使用 dtype=np.float32 转换数组。如果您不能使用简单精度或混合精度，另一种昂贵的解决方案可能是使用专用于此的专业 GPU。
此外，多个表达式可以预先计算并存储在常量中。例如，这包括 array.astype(np.float32)、math.cos(angi) 和 math.sin(angi)。其他一些表达式可以存储在临时变量中，因为编译器可能无法有效地分解代码（主要是因为 trigonometric functions）。
此外，三角函数通常非常昂贵，尤其是当您希望计算符合 IEEE-754 标准时（1.0/SDD 调用很可能就是这种情况）。您可以改用近似值。 CUDA 提供了 math.xxx、__cosf 和 __sinf 内在函数，它们应该更快（但如果使用它们，请注意结果）。我不确定您是否可以直接调用它们，但您可以将参数 __tanf 添加到 JIT 装饰器中，它可以为您执行此操作。
我认为使用 32x8 的 2D 线程块可能会快一些，因为线程被打包在包含 32 个线程和 GPU 的 warp 中。但最好的解决方案是检查许多不同块大小的性能。
如果所有这些还不够，您可以尝试使用共享内存来减少每个块执行的指令量，因为每个块对某些表达式进行多次重新计算。
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
