1

I would like to implement the code below avoiding for loops to increase speed. Is there any way to do this so that I can create a sphere centered in the numpy array?

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

Adopting the solution from the recommended link Python vectorizing nested for loops I was able to replace the code. I am getting unexplained artifacts in the final plot however. There are rings appearing around the central sphere. What could be causing these?

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

This code give a sphere centered with some ring artifacts around Matplotlib imshow plot

john
  • 57
  • 1
  • 7
  • 2
    This should be very close - https://stackoverflow.com/questions/39667089/python-vectorizing-nested-for-loops – Divakar Oct 05 '20 at 16:57

1 Answers1

0

I realized that using unsigned int was causing errors in subtraction. The final working solution is below

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
john
  • 57
  • 1
  • 7