Created
December 22, 2015 21:14
-
-
Save adler-j/94ea8ee27a6a11d58a14 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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