yade-users team mailing list archive

[Karol Brzezinski <question696923@xxxxxxxxxxxxxxxxxxxxx>]

Question #696923 on Yade changed:
https://answers.launchpad.net/yade/+question/696923

Karol Brzezinski proposed the following answer:
Hi,

please find my script below. It takes about 7 minutes on my computer to
take 10 slices with 0.6 Mpx resolution (pixel size 0.000005 m), but
obviously you need more slices for your reconstruction.

Cheers,
KB

#######################
import numpy as np
from PIL import Image

n_slices = 50

x_s = np.array([])
y_s = np.array([])
z_s = np.array([])
radii = np.array([])
for b in O.bodies:
	if isinstance(b.shape, Sphere):
		pos = b.state.pos
		x_s = np.append(x_s,pos[0])
		y_s = np.append(y_s,pos[1])
		z_s = np.append(z_s,pos[2])
		radii = np.append(radii,b.shape.radius)
		
px_size = 0.000005 # pixel size (m)
(x_min,y_min,z_min),(x_max,y_max,z_max) = aabbExtrema()
real_h = z_max-z_min
real_w = x_max - x_min # I assume here that slices will be perpendicular to y-axis

img_h = ceil(real_h/px_size) # image height in px
img_w = ceil(real_w/px_size) # image height in px

"""
The idea is as follows:
- divide the cross section into squares (representing circles)
- for the center of every square (pixel) check whether any center of the sphere is closer than its radius
- if True color the pixel white, otherwise black

One could probabily optimize this. I only took adance on the vectorisation in numpy. 
I have made a "manual" decision for this algorithm: since the number of pixels is greater than number of bodies I loop over spheres and do vectorized operation on the "pixel" array.

"""

x_px = np.tile(np.arange(x_min+px_size/2,x_min+px_size*img_w,px_size),img_h)
z_px = np.repeat(np.arange(z_min+px_size/2,z_min+px_size*img_h,px_size),img_w)

for y_slice in np.linspace(y_min,y_max,n_slices+2)[1:-1]:# take n_slices slices (linspace prepares n_slices+2, but I ommit first and last since it is not interesting)
	y_px = np.repeat(y_slice,img_h*img_w)#
	pixels_state = np.zeros(img_h*img_w).astype(bool)

	count = 0
	for i in range(len(radii)):
		radius = radii[i]
		x = x_s[i]
		y = y_s[i]
		z = z_s[i]
		if abs(y-y_slice)<=radius:# I make use of simplified assumption, I only check pixels for the spheres that are close enough to my plane.
			pixels_state_tmp = np.sqrt((x-x_px)**2+(y-y_px)**2+(z-z_px)**2)<radius
			pixels_state = pixels_state + pixels_state_tmp
		
	img_array = 256*pixels_state.reshape(img_h,img_w).astype(np.uint8)
			
	img = Image.fromarray(img_array)	
	img.save('slice_y_{:.5f}.gif'.format(y_slice))

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.