在没有使用网格网格的 for 循环的情况下在数组中心创建球体会创建外壳工件

Question

我想实现下面的代码，避免 for 循环以提高速度。有什么办法可以创建一个以 numpy 数组为中心的球体？

def Create_Sphere(square_numpy_array, pixel_min, pixel_max, HU, Radius_sq ):


new_array = np.empty_like(square_numpy_array)

for k in range(pixel_min, pixel_max, 1):
    for i in range(pixel_min, pixel_max, 1):
        for j in range(pixel_min, pixel_max, 1):
            r_sq = (i - 255)**2 + (j - 255)**2 + (k - 255)**2
            if r_sq <= Radius_sq:
                new_array[k, i, j] = HU + 1000
return new_array

采用推荐链接Python vectorizing nested for loops 中的解决方案，我能够替换代码。然而，我在最终情节中得到了无法解释的伪影。中心球体周围出现环。这些可能是什么原因造成的？

def Create_Sphere_CT(HU=12):

     radius = np.uint16(100) #mm
     Radius_sq_pixels = np.uint16((radius *2)**2 )
     sphere_pixel_HU = np.uint16(HU + 1000) #dtype controlled for memory
     center_pixel = np.uint16(400/2-1)  
     new_array = np.zeros((400,400,400), dtype = np.uint16)

     m,n,r = new_array.shape
     x = np.arange(0, m, 1, dtype = np.uint16)
     y = np.arange(0, n, 1, dtype = np.uint16)
     z = np.arange(0, r, 1, dtype = np.uint16)

     xx,yy,zz = np.meshgrid(x,y,z, indexing = 'ij')

     X = (xx - center_pixel)
     xx = None #free memory once variable is used
     Y = (yy - center_pixel)
     yy= None #free memory once variable is used
     Z = (zz - center_pixel)
     zz = None#free memory once variable is used

     mask = (X**2 + Y**2 + Z**2) < Radius_sq_pixels  #create sphere mask
     new_array = sphere_pixel_HU * mask  #assign values 

     return new_array

这段代码给出了一个以一些环形伪影为中心的球体

Answer 1

我意识到使用 unsigned int 会导致减法错误。最终的工作解决方案如下

def Sphere(HU):
num_pix = int(400)
radius =100 
Radius_sq_pixels = (radius)**2
sphere_pixel_HU = HU 
center_pixel = int(num_pix/2-1) 
new_array = np.zeros((num_pix, num_pix, num_pix))

m,n,r = new_array.shape
x = np.arange(0, m, 1)
y = np.arange(0, n, 1)
z = np.arange(0, r, 1)

xx,yy,zz = np.meshgrid(x,y,z,indexing = 'ij',sparse = True)
X = (xx - center_pixel)
Y = (yy - center_pixel)
Z = (zz - center_pixel)

mask = ((X**2) + (Y**2) + (Z**2)) < Radius_sq_pixels  #create sphere mask

new_array = sphere_pixel_HU * mask  #assign values 
new_array = new_array.astype(np.uint16) #change datatype

plt.imshow(new_array[int(num_pix/2)])
plt.show()

return new_array