Skip to content

Instantly share code, notes, and snippets.

@adler-j
Created December 22, 2015 21:14
Show Gist options
  • Save adler-j/94ea8ee27a6a11d58a14 to your computer and use it in GitHub Desktop.
Save adler-j/94ea8ee27a6a11d58a14 to your computer and use it in GitHub Desktop.
277 # Blank image
278 1 71 71.0 0.0 p = np.zeros(space.shape)
279
280 # Create the pixel grid
281 1 308 308.0 0.0 grid_in = space.grid.meshgr
id()
282 1 25 25.0 0.0 minp = space.grid.min()
283 1 18 18.0 0.0 maxp = space.grid.max()
284
285 # move points to [-1, 1]
286 1 4 4.0 0.0 grid = []
287 4 21 5.2 0.0 for i in range(3):
288 3 23 7.7 0.0 meani = (minp[i] + maxp
[i]) / 2.0
289 3 18 6.0 0.0 diffi = (maxp[i] - minp
[i]) / 2.0
290 3 93 31.0 0.0 grid += [(grid_in[i] -
meani) / diffi]
291
292 11 126 11.5 0.0 for ellip in ellipses:
293 10 48 4.8 0.0 I = ellip[0]
294 10 99 9.9 0.0 a2 = ellip[1] ** 2
295 10 54 5.4 0.0 b2 = ellip[2] ** 2
296 10 53 5.3 0.0 c2 = ellip[3] ** 2
297 10 46 4.6 0.0 x0 = ellip[4]
298 10 41 4.1 0.0 y0 = ellip[5]
299 10 40 4.0 0.0 z0 = ellip[6]
300 10 122 12.2 0.0 phi = ellip[7] * np.pi
/ 180
301 10 56 5.6 0.0 theta = ellip[8] * np.p
i / 180
302 10 53 5.3 0.0 psi = ellip[9] * np.pi
/ 180
303
304 10 108 10.8 0.0 scales = [1 / a2, 1 / b
2, 1 / c2]
305
306 # Create the offset x,y
and z values for the grid
307 10 97 9.7 0.0 if any([phi, theta, psi
]):
308 # Calculate the poi
nts that could possibly be inside the volume
309 # Since the points
are rotated, we cannot do anything directional
310 # without more logi
c
311 2 196 98.0 0.0 center = (np.array(
[x0, y0, z0]) + 1.0) / 2.0
312 2 376 188.0 0.0 max_radius = np.sqr
t(np.max([a2, b2, c2]))
313 2 466 233.0 0.0 idx, shapes = _gets
hapes(center, max_radius, space.shape)
314
315 # Rotate the points
to the expected coordinate system.
316 2 33 16.5 0.0 cphi = np.cos(phi)
317 2 24 12.0 0.0 sphi = np.sin(phi)
318 2 17 8.5 0.0 ctheta = np.cos(the
ta)
319 2 16 8.0 0.0 stheta = np.sin(the
ta)
320 2 17 8.5 0.0 cpsi = np.cos(psi)
321 2 16 8.0 0.0 spsi = np.sin(psi)
322
323 2 17 8.5 0.0 mat = np.array([[cp
si * cphi - ctheta * sphi * spsi,
324 2 14 7.0 0.0 cp
si * sphi + ctheta * cphi * spsi,
325 2 10 5.0 0.0 sp
si * stheta],
326 2 15 7.5 0.0 [-s
psi * cphi - ctheta * sphi * cpsi,
327 2 12 6.0 0.0 -s
psi * sphi + ctheta * cphi * cpsi,
328 2 8 4.0 0.0 cp
si * stheta],
329 2 8 4.0 0.0 [st
heta * sphi,
330 2 9 4.5 0.0 -s
theta * cphi,
331 2 54 27.0 0.0 ct
heta]])
332
333 8 73 9.1 0.0 subgrid = [g[idi] f
or g, idi in zip(grid, shapes)]
334 2 8 4.0 0.0 offset_points = [ve
c * (xi - x0i)[..., np.newaxis]
335 2 10 5.0 0.0 fo
r xi, vec, x0i in zip(subgrid,
336 2 16 8.0 0.0
mat.T,
337 8 373 46.6 0.0
[x0, y0, z0])]
338 2 158754 79377.0 11.9 rotated = offset_po
ints[0] + offset_points[1] + offset_points[2]
339 2 24685 12342.5 1.9 np.square(rotated,
out=rotated)
340 2 346955 173477.5 26.1 radius = np.dot(rot
ated, scales)
341 else:
342 # Calculate the poi
nts that could possibly be inside the volume
343 8 442 55.2 0.0 center = (np.array(
[x0, y0, z0]) + 1.0) / 2.0
344 8 211 26.4 0.0 max_radius = np.sqr
t([a2, b2, c2])
345 8 1310 163.8 0.1 idx, shapes = _gets
hapes(center, max_radius, space.shape)
346
347 32 253 7.9 0.0 subgrid = [g[idi] f
or g, idi in zip(grid, shapes)]
348 8 39 4.9 0.0 squared_dist = [ai
* (xi - x0i)**2
349 8 35 4.4 0.0 for
xi, ai, x0i in zip(subgrid,
350 8 32 4.0 0.0
scales,
351 32 870 27.2 0.1
[x0, y0, z0])]
352
353 8 254687 31835.9 19.1 radius = squared_di
st[0] + squared_dist[1] + squared_dist[2]
354
355 # Find the pixels withi
n the ellipse
356 10 91403 9140.3 6.9 inside = radius <= 1
357
358 # Add the ellipse inten
sity to those pixels
359 10 416653 41665.3 31.3 p[idx][inside] += I
360
361 1 32138 32138.0 2.4 return space.element(p)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment