forked from abnowack/pytracer
-
Notifications
You must be signed in to change notification settings - Fork 1
/
raytrace.py
executable file
·108 lines (82 loc) · 3.45 KB
/
raytrace.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
from mesh import create_hollow, create_rectangle, create_circle
from material import Material
from solid import Solid
from simulation import Simulation
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.interpolation import rotate
import sys
# TODO : Normalization isn't correct
def inverse_radon(radon, thetas):
"""
Reconstruct using Filtered Back Projection.
Weighting assumes thetas are equally spaced
radon size must be even
"""
pad_value = int(2 ** (np.ceil(np.log(2 * np.size(radon, 0)) / np.log(2))))
pre_pad = int((pad_value - len(radon[:, 0])) / 2)
post_pad = pad_value - len(radon[:, 0]) - pre_pad
f = np.fft.fftfreq(pad_value)
ramp_filter = 2. * np.abs(f)
reconstruction_image = np.zeros((np.size(radon, 0), np.size(radon, 0)))
for i, theta in enumerate(thetas):
filtered = np.real(np.fft.ifft(np.fft.fft(np.pad(radon[:, i], (pre_pad, post_pad), 'constant', constant_values=(0, 0))) * ramp_filter))[pre_pad:-post_pad]
back_projection = rotate(np.tile(filtered, (np.size(radon, 0), 1)), theta, reshape=False, mode='constant')
reconstruction_image += back_projection * 2 * np.pi / len(thetas)
return reconstruction_image
def build_shielded_geometry():
air = Material(0.1, color='white')
u235_metal = Material(1.0, color='green')
poly = Material(1.0, color='red')
steel = Material(1.0, color='orange')
box = create_hollow(create_rectangle(20., 10.), create_rectangle(18., 8.))
hollow_circle = create_hollow(create_circle(3.9), create_circle(2.9))
hollow_circle.translate([-9 + 3.9 + 0.1, 0.])
small_box_1 = create_rectangle(2., 2.)
small_box_1.translate([6., 2.])
small_box_2 = create_rectangle(2., 2.)
small_box_2.translate([6., -2.])
#sim = Simulation(air, 50., 45., 'arc')
sim = Simulation(air, 100, diameter=50., detector='plane', detector_width=30.)
sim.detector.width = 30.
sim.geometry.solids.append(Solid(box, steel, air))
sim.geometry.solids.append(Solid(hollow_circle, steel, air))
sim.geometry.solids.append(Solid(small_box_1, poly, air))
sim.geometry.solids.append(Solid(small_box_2, steel, air))
sim.geometry.flatten()
return sim
def ray_trace_test_geometry():
air = Material(0.0, color='white')
steel = Material(1.0, color='red')
box = create_hollow(create_rectangle(12., 12.), create_rectangle(10., 10.))
ring = create_hollow(create_circle(12.), create_circle(10.))
box.rotate(45.)
sim = Simulation(air, diameter=50.)
sim.detector.width = 30.
sim.geometry.solids.append(Solid(ring, steel, air))
sim.geometry.flatten()
return sim
def main():
sim = build_shielded_geometry()
plt.figure()
sim.draw(True)
n_angles = 100
angles = np.linspace(0., 180., n_angles + 1)[:-1]
radon = sim.radon_transform(angles)
plt.figure()
plt.imshow(radon, cmap=plt.cm.Greys_r, interpolation='none',
aspect='auto')
plt.xlabel('Angle')
plt.ylabel('Radon Projection')
plt.colorbar()
plt.figure()
recon_image = inverse_radon(radon, angles)
extent = [-sim.detector.width / 2., sim.detector.width / 2.]
plt.imshow(recon_image.T[:, ::-1], interpolation='none', extent=extent *
2)
plt.xlabel('X (cm)')
plt.ylabel('Y (cm)')
plt.colorbar()
plt.show()
if __name__ == "__main__":
sys.exit(int(main() or 0))