WIP

Author: Georg Schramm

README.md

@@ -1,41 +1,35 @@
-# PETSIRD v0.2 - parallelproj repository
-
-The purpose of this repo is to provide a starting point for developing software that uses PETSIRD.
-Here we will focus on python usecases using parallelproj projection library.
+# Analytic simulation of petsird v0.2 data using parallelproj
 
 ## Create your conda / mamba environment
 
 ```
 conda env create -f environment.yml
-```
-
-## Run a first PETSIRD example
-
-```
 conda activate petsird-v0.2-parallelproj
 cd python
-python petsird_generator.py > test.bin
-python petsird_analysis.py < test.bin
 ```
 
-## Projects to work on
-
-### Read scanner geometry from a PETSIRD file and visualize it
+## Simulate petsird LM data
 
-**TODO**: 
-- start from this [parallel example here](https://parallelproj.readthedocs.io/en/stable/auto_examples/01_pet_geometry/03_run_block_scanner.html#sphx-glr-auto-examples-01-pet-geometry-03-run-block-scanner-py) and this [petsird example here](python/petsird_plot_scanner_and_events.py)
-- map PETSIRD scanner detector elements to LOR endpoints
-- How to map BoxShape to a single LOR point? (center of the Box is not want we want. We need most likely interaction point.)
-- use LOR endpoints to setup [BlockScannerModules](https://parallelproj.readthedocs.io/en/stable/python_api.html#parallelproj.pet_scanners.BlockPETScannerModule) tha can be used in [ModularizedPETScannerGeometry](https://parallelproj.readthedocs.io/en/stable/python_api.html#parallelproj.pet_scanners.ModularizedPETScannerGeometry) or a [RegularPolygonPETScannerGeometry](https://parallelproj.readthedocs.io/en/stable/python_api.html#parallelproj.pet_scanners.RegularPolygonPETScannerGeometry) - depending of the system geometry and module definition. Using the latter is more convient for the definition of sinograms, the former is better for systems with less symmetries where it is less straight forward to define "sinograms".
-
-### Simulation of PETSIRD LM data using parallelproj
+```
+python 01_analytic_petsird_lm_simulator.py
+```
 
-Start from this [parallelproj example](https://parallelproj.readthedocs.io/en/stable/auto_examples/06_listmode_algorithms/01_listmode_mlem.html#sphx-glr-auto-examples-06-listmode-algorithms-01-listmode-mlem-py).
+The simulation can be customized in many ways (number of counts to simulate,
+uniform or non-uniform efficiencies ...) via command line options.
 
-**TODO**: 
-- rewrite (simplify) util function that converts sinogram to listmode including event detectors IDs instead of coordinates. Write an updated noisy emission sinogram to LM converter similar to  [this](https://parallelproj.readthedocs.io/en/stable/_modules/parallelproj/projectors.html#RegularPolygonPETProjector.convert_sinogram_to_listmode) but skipping the detector ID to world coordinate lookup.
-- write scanner geometry to file
+These option can be listed via
+```
+python 01_analytic_petsird_lm_simulator.py
+```
 
+## Run a listmode OSEM recon on the simulated
 
+```
+python 02_recon_block_scanner_listmode_data.py
+```
 
+Thes command line optione for the LM OSEM recon script can be listed via
+```
+python 02_recon_block_scanner_listmode_data.py -h
+```
 

python/01_analytic_petsird_lm_simulator.py

@@ -44,7 +44,6 @@ def module_pair_eff_from_sgd(i_sgd: int, uniform: bool = False) -> float:
     return res
 
 
-# parse the command line
 def parse_int_tuple(arg):
     return tuple(map(int, arg.split(",")))
 
@@ -63,7 +62,7 @@ def parse_float_tuple(arg):
 
 parser.add_argument("--fname", type=str, default="simulated_lm_file.bin")
 parser.add_argument("--output_dir", type=str, default="my_lm_sim")
-parser.add_argument("--num_true_counts", type=int, default=int(1e6))
+parser.add_argument("--num_true_counts", type=int, default=int(4e6))
 parser.add_argument("--skip_plots", action="store_true")
 parser.add_argument("--check_backprojection", default=False, action="store_true")
 parser.add_argument("--num_epochs_mlem", type=int, default=0)
@@ -78,7 +77,7 @@ def parse_float_tuple(arg):
 parser.add_argument(
     "--phantom",
     type=str,
-    default="uniform_cylinder",
+    default="squares",
     choices=["uniform_cylinder", "squares"],
 )
 

python/recon_block_scanner_listmode_data.py renamed to python/02_recon_block_scanner_listmode_data.py

@@ -76,10 +76,18 @@ def transform_BoxShape(
     )
 
 
+def parse_int_tuple(arg):
+    return tuple(map(int, arg.split(",")))
+
+
+def parse_float_tuple(arg):
+    return tuple(map(float, arg.split(",")))
+
+
 # %%
 
 parser = argparse.ArgumentParser()
-parser.add_argument("--fname", type=str, default="sim_lm.bin")
+parser.add_argument("--fname", type=str, default="my_lm_sim/simulated_lm_file.bin")
 parser.add_argument("--num_epochs", type=int, default=5)
 parser.add_argument("--num_subsets", type=int, default=20)
 parser.add_argument("--img_shape", type=parse_int_tuple, default=(100, 100, 11))
@@ -223,13 +231,13 @@ def transform_BoxShape(
                 proj.event_tofbins = xp.full(
                     start_coords.shape[0], tofbin, dtype="int32"
                 )
+                proj.tof = True
                 sens_img += proj.adjoint(start_el_eff * end_el_eff * module_pair_eff)
 
 print("")
 
 # for some reason we have to divide the sens image by the number of TOF bins
 # right now unclear why that is
-sens_img /= num_tofbins
 sens_img = res_model.adjoint(sens_img)
 
 # %%
@@ -294,6 +302,7 @@ def transform_BoxShape(
     )
     lm_subset_projs[i_subset].tof_parameters = tof_params
     lm_subset_projs[i_subset].event_tofbins = tof_bin[sl]
+    lm_subset_projs[i_subset].tof = True
 
 for i_epoch in range(num_epochs):
     for i_subset, sl in enumerate(subset_slices):