add noise + event generation and plot of first 3 events

python/parallelproj_sim.py

@@ -9,7 +9,7 @@
 # are performed on the cpu or cuda gpu
 
 dev = "cpu"
-expected_num_trues = 1e6
+expected_num_trues = 5e6
 np.random.seed(1)
 
 # ----------------------------------------------------------------------------
@@ -63,51 +63,68 @@
 # calculate the sensitivity image
 sens_img = projector.adjoint(np.ones(noise_free_sinogram.shape, device=dev, dtype=np.float32))
 
-## get the start and end points of all sinogram bins
-#xstart, xend = lor_descriptor.get_lor_coordinates()
+# get a lookup table that contains the world coordinates of all scanner detectors
+scanner_lut = lor_descriptor.scanner.all_lor_endpoints
 
 # get the two dimensional indices of all sinogram bins
 start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices()
 
-# convert the 2D indices to a 1D index
-det_start_index = lor_descriptor.scanner.num_modules * start_mods + start_inds
-det_end_index = lor_descriptor.scanner.num_modules * end_mods + end_inds
+# generate two sinograms that contain the linearized detector start and end indices
+sino_det_start_index =  lor_descriptor.scanner.num_lor_endpoints_per_module[0] * start_mods + start_inds
+sino_det_end_index =  lor_descriptor.scanner.num_lor_endpoints_per_module[0] * end_mods + end_inds
 
 # add poisson noise to the noise free sinogram
 noisy_sinogram = np.random.poisson(noise_free_sinogram)
 
 # ravel the noisy sinogram and the detector start and end "index" sinograms 
 noisy_sinogram = noisy_sinogram.ravel()
-det_start_index = det_start_index.ravel()
-det_end_index = det_end_index.ravel()
+sino_det_start_index = sino_det_start_index.ravel()
+sino_det_end_index = sino_det_end_index.ravel()
+
+# ----------------------------------------------------------------------------
+# ----------------------------------------------------------------------------
 
 # generate listmode data from the noisy sinogram
-inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram)
+event_sino_inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram)
+# shuffle the event sinogram indices
+np.random.shuffle(event_sino_inds)
 
+event_det_id_1 = sino_det_start_index[event_sino_inds]
+event_det_id_2 = sino_det_end_index[event_sino_inds]
 
 # ----------------------------------------------------------------------------
 # ----------------------------------------------------------------------------
 # ----------------------------------------------------------------------------
 # ----------------------------------------------------------------------------
 
 # show the scanner geometry and one view in one sinogram plane
-fig = plt.figure()
+fig = plt.figure(figsize=(8, 8))
 ax = fig.add_subplot(projection="3d")
-lor_descriptor.scanner.show_lor_endpoints(ax)
-lor_descriptor.show_views(
-    ax,
-    views=np.asarray([lor_descriptor.num_views // 4], device=dev),
-    planes=np.asarray([lor_descriptor.scanner.num_modules // 2], device=dev),
-)
+lor_descriptor.scanner.show_lor_endpoints(ax, show_linear_index=True, annotation_fontsize=4)
+
+# plot the LORs of the first n events
+for i, col in enumerate(('red','green','blue')):
+    xstart = scanner_lut[event_det_id_1[i], :]
+    xend = scanner_lut[event_det_id_2[i], :]
+
+    ax.plot(
+        [xstart[0], xend[0]],
+        [xstart[1], xend[1]],
+        [xstart[2], xend[2]],
+        color=col,
+        linewidth=1.,
+    )
+
 fig.tight_layout()
 fig.show()
-
-fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5))
-ax2[0].imshow(np.asarray(to_device(img[:, :, 3], "cpu")))
-if projector.tof:
-    ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu")))
-else:
-    ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu")))
-ax2[2].imshow(np.asarray(to_device(sens_img[:, :, 3], "cpu")))
-fig2.tight_layout()
-fig2.show()
+#
+#fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5))
+#ax2[0].imshow(np.asarray(to_device(img[:, :, 3], "cpu")))
+#if projector.tof:
+#    ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu")))
+#else:
+#    ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu")))
+#ax2[2].imshow(np.asarray(to_device(sens_img[:, :, 3], "cpu")))
+#fig2.tight_layout()
+#fig2.show()
+#