diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index e775cd1..ea2e1b5 100644 --- a/python/parallelproj_sim.py +++ b/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,27 +63,34 @@ # 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] # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -91,23 +98,33 @@ # ---------------------------------------------------------------------------- # 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() +# \ No newline at end of file