add script to produce attenuation image

Author: Georg Schramm

python/01_analytic_petsird_lm_simulator.py

@@ -70,55 +70,127 @@ def parse_float_tuple(arg):
 # %%
 # parse the command line for the input parameters below
 parser = argparse.ArgumentParser(
-    description="Analytic simulation of PETSIRD v0.7.2 listmode data for a block PET scanner"
+    description="Analytic simulation of PETSIRD v0.7.2 listmode data for a block PET scanner",
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 )
 
 # Output configuration group
-output_group = parser.add_argument_group('Output Configuration')
-output_group.add_argument("--fname", type=str, default="simulated_petsird_lm_file.bin", help="name of the output LM file")
-output_group.add_argument("--output_dir", type=str, default=None, help="directory to save output files")
-output_group.add_argument("--skip_writing", default=False, action="store_true", help="skip writing the LM data to a file")
+output_group = parser.add_argument_group("Output Configuration")
+output_group.add_argument(
+    "--fname",
+    type=str,
+    default="simulated_petsird_lm_file.bin",
+    help="name of the output LM file",
+)
+output_group.add_argument(
+    "--output_dir", type=str, default=None, help="directory to save output files"
+)
+output_group.add_argument(
+    "--skip_writing",
+    default=False,
+    action="store_true",
+    help="skip writing the LM data to a file",
+)
 
 # Image and phantom configuration group
-image_group = parser.add_argument_group('Image and Phantom Configuration')
+image_group = parser.add_argument_group("Image and Phantom Configuration")
 image_group.add_argument(
     "--phantom",
     type=str,
     default="cylinder",
     choices=["cylinder", "uniform_cylinder", "squares", "points"],
     help="phantom to simulate",
 )
-image_group.add_argument("--img_shape", type=parse_int_tuple, default=(55, 55, 19), help="shape of the image to simulate")
-image_group.add_argument("--voxel_size", type=parse_float_tuple, default=(2.0, 2.0, 2.0), help="voxel size in mm used in the simulation")
+image_group.add_argument(
+    "--img_shape",
+    type=parse_int_tuple,
+    default=(55, 55, 19),
+    help="shape of the image to simulate",
+)
+image_group.add_argument(
+    "--voxel_size",
+    type=parse_float_tuple,
+    default=(2.0, 2.0, 2.0),
+    help="voxel size in mm used in the simulation",
+)
 
 # Physics and resolution parameters group
-physics_group = parser.add_argument_group('Physics / simulation parameters')
-physics_group.add_argument("--fwhm_mm", type=float, default=2.5, help="FWHM of the image space resolution model in mm")
-physics_group.add_argument("--tof_fwhm_mm", type=float, default=20.0, help="FWHM of the TOF resolution model in mm")
-physics_group.add_argument("--num_true_counts", type=int, default=int(4e6), help="number of true coincidences to simulate")
+physics_group = parser.add_argument_group("Physics / simulation parameters")
+physics_group.add_argument(
+    "--fwhm_mm",
+    type=float,
+    default=2.5,
+    help="FWHM of the image space resolution model in mm",
+)
+physics_group.add_argument(
+    "--tof_fwhm_mm",
+    type=float,
+    default=20.0,
+    help="FWHM of the TOF resolution model in mm",
+)
+physics_group.add_argument(
+    "--num_true_counts",
+    type=int,
+    default=int(4e6),
+    help="number of true coincidences to simulate",
+)
 
 # Efficiency parameters group
-efficiency_group = parser.add_argument_group('Detection Efficiency Parameters')
-efficiency_group.add_argument("--uniform_crystal_eff", action="store_true", help="use uniform crystal efficiencies, otherwise use a random distribution")
-efficiency_group.add_argument("--uniform_sg_eff", action="store_true", help="use uniform symmetry group (module pair) efficiencies, otherwise pseudo random pattern is used")
+efficiency_group = parser.add_argument_group("Detection Efficiency Parameters")
+efficiency_group.add_argument(
+    "--uniform_crystal_eff",
+    action="store_true",
+    help="use uniform crystal efficiencies, otherwise use a random distribution",
+)
+efficiency_group.add_argument(
+    "--uniform_sg_eff",
+    action="store_true",
+    help="use uniform symmetry group (module pair) efficiencies, otherwise pseudo random pattern is used",
+)
 
 # Time block configuration group
-time_group = parser.add_argument_group('Time Block Configuration')
-time_group.add_argument("--num_time_blocks", type=int, default=3, help="number of time blocks to split the data into")
-time_group.add_argument("--event_block_duration", type=int, default=100, help="duration of each time block in ms")
+time_group = parser.add_argument_group("Time Block Configuration")
+time_group.add_argument(
+    "--num_time_blocks",
+    type=int,
+    default=3,
+    help="number of time blocks to split the data into",
+)
+time_group.add_argument(
+    "--event_block_duration",
+    type=int,
+    default=100,
+    help="duration of each time block in ms",
+)
 
 # Reconstruction and analysis group
-recon_group = parser.add_argument_group('Reconstruction and Analysis')
-recon_group.add_argument("--num_epochs_mlem", type=int, default=0, help="number of epochs for MLEM reconstruction of histogrammed data, 0 means no MLEM reconstruction")
-recon_group.add_argument("--check_backprojection", default=False, action="store_true", help="check the backprojection of the TOF histogram and LM data")
+recon_group = parser.add_argument_group("Reconstruction and Analysis")
+recon_group.add_argument(
+    "--num_epochs_mlem",
+    type=int,
+    default=0,
+    help="number of epochs for MLEM reconstruction of histogrammed data, 0 means no MLEM reconstruction",
+)
+recon_group.add_argument(
+    "--check_backprojection",
+    default=False,
+    action="store_true",
+    help="check the backprojection of the TOF histogram and LM data",
+)
 
 # Visualization and debugging group
-visual_group = parser.add_argument_group('Visualization and Debugging')
-visual_group.add_argument("--skip_plots", action="store_true", help="skip plotting the scanner geometry and TOF profile")
+visual_group = parser.add_argument_group("Visualization and Debugging")
+visual_group.add_argument(
+    "--skip_plots",
+    action="store_true",
+    help="skip plotting the scanner geometry and TOF profile",
+)
 
 # General options group
-general_group = parser.add_argument_group('General Options')
-general_group.add_argument("--seed", type=int, default=0, help="random seed for reproducibility")
+general_group = parser.add_argument_group("General Options")
+general_group.add_argument(
+    "--seed", type=int, default=0, help="random seed for reproducibility"
+)
 
 args = parser.parse_args()
 

python/02_reconstruct_petsird.py

@@ -3,6 +3,7 @@
 from pathlib import Path
 import parallelproj
 import petsird
+import hashlib
 
 # for the parallelproj recon we can use cupy as array backend if available
 # otherwise we fall back to numpy
@@ -21,13 +22,44 @@
 )
 
 
+def get_params_hash(
+    img_shape: tuple[int, int, int],
+    voxel_size: tuple[float, float, float],
+    tof: bool,
+    att_img: np.ndarray | None = None,
+) -> str:
+    """Generate a hash string for all sensitivity image parameters."""
+    # Create a hash based on all parameters that affect sensitivity image calculation
+    hash_obj = hashlib.sha256()
+
+    # Add image shape
+    hash_obj.update(str(img_shape).encode())
+
+    # Add voxel size (convert to string with fixed precision to avoid floating point issues)
+    voxel_size_str = f"{voxel_size[0]:.6f},{voxel_size[1]:.6f},{voxel_size[2]:.6f}"
+    hash_obj.update(voxel_size_str.encode())
+
+    # Add TOF flag
+    hash_obj.update(str(tof).encode())
+
+    # Add attenuation image data if present
+    if att_img is not None:
+        hash_obj.update(att_img.tobytes())
+    else:
+        hash_obj.update(b"no_attenuation")
+
+    # Return first 12 characters to reduce collision probability with more parameters
+    return hash_obj.hexdigest()[:6]
+
+
 # %%
 ################################################################################
 #### PARSE THE COMMAND LINE ####################################################
 ################################################################################
 
 parser = argparse.ArgumentParser(
-    description="PETSIRD analytic simulator reconstruction"
+    description="PETSIRD analytic simulator reconstruction",
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 )
 parser.add_argument("fname", type=str, help="Path to the PETSIRD listmode file")
 parser.add_argument(
@@ -210,7 +242,10 @@
     print("Using ones as sensitivity image ...")
     sens_img = xp.ones(img_shape, dtype="float32")
 else:
-    sens_img_path = Path(fname).with_suffix(".sens_img.npy")
+    # Generate hash for all sensitivity parameters
+    params_hash = get_params_hash(img_shape, voxel_size, tof, att_img)
+
+    sens_img_path = Path(fname).with_name(f"{Path(fname).stem}_sens_{params_hash}.npy")
 
     if sens_img_path.exists():
         print(f"Loading sensitivity image from {sens_img_path}")
@@ -235,25 +270,25 @@
 
     xp.save(sens_img_path, sens_img)
 
-# %%
-################################################################################
-### CHECK WHETHER THE CALC SENS IMAGE MATCHES THE REF. SENSE IMAGE #############
-################################################################################
-
-ref_sens_img_path = Path(fname).parent / "reference_sensitivity_image.npy"
-
-if ref_sens_img_path.exists():
-    print(f"loading reference sensitivity image from {ref_sens_img_path}")
-    ref_sens_img = np.load(ref_sens_img_path)
-    if ref_sens_img.shape == sens_img.shape:
-        if np.allclose(np.asarray(sens_img), ref_sens_img):
-            print(
-                f"calculated sensitivity image matches reference image {ref_sens_img_path}"
-            )
-        else:
-            print(
-                f"calculated sensitivity image does NOT match reference image {ref_sens_img_path}"
-            )
+## %%
+#################################################################################
+#### CHECK WHETHER THE CALC SENS IMAGE MATCHES THE REF. SENSE IMAGE #############
+#################################################################################
+#
+# ref_sens_img_path = Path(fname).parent / "reference_sensitivity_image.npy"
+#
+# if ref_sens_img_path.exists():
+#    print(f"loading reference sensitivity image from {ref_sens_img_path}")
+#    ref_sens_img = np.load(ref_sens_img_path)
+#    if ref_sens_img.shape == sens_img.shape:
+#        if np.allclose(np.asarray(sens_img), ref_sens_img):
+#            print(
+#                f"calculated sensitivity image matches reference image {ref_sens_img_path}"
+#            )
+#        else:
+#            print(
+#                f"calculated sensitivity image does NOT match reference image {ref_sens_img_path}"
+#            )
 
 # %%
 ################################################################################
@@ -325,25 +360,32 @@
 
 non_tof_backproj = proj.adjoint(xp.ones(coords0.shape[0], dtype="float32"))
 
-xp.save(Path(fname).with_suffix(".non_tof_backproj.npy"), non_tof_backproj)
+non_tof_backproj_path = Path(fname).with_name(
+    f"{Path(fname).stem}_nontof_backproj_{params_hash}.npy"
+)
+xp.save(non_tof_backproj_path, non_tof_backproj)
 
 #### HACK assumes same TOF parameters for all module type pairs
-sigma_tof = scanner_info.tof_resolution[0][0] / 2.35
-tof_bin_edges = scanner_info.tof_bin_edges[0][0].edges
-num_tofbins = tof_bin_edges.size - 1
-tofbin_width = float(tof_bin_edges[1] - tof_bin_edges[0])
-
-tof_params = parallelproj.TOFParameters(
-    num_tofbins=num_tofbins, tofbin_width=tofbin_width, sigma_tof=sigma_tof
-)
+if tof:
+    sigma_tof = scanner_info.tof_resolution[0][0] / 2.35
+    tof_bin_edges = scanner_info.tof_bin_edges[0][0].edges
+    num_tofbins = tof_bin_edges.size - 1
+    tofbin_width = float(tof_bin_edges[1] - tof_bin_edges[0])
+
+    tof_params = parallelproj.TOFParameters(
+        num_tofbins=num_tofbins, tofbin_width=tofbin_width, sigma_tof=sigma_tof
+    )
 
-proj.tof_parameters = tof_params
-proj.event_tofbins = xp.asarray(signed_tof_bins).copy()
-proj.tof = True
+    proj.tof_parameters = tof_params
+    proj.event_tofbins = xp.asarray(signed_tof_bins).copy()
+    proj.tof = True
 
-tof_backproj = proj.adjoint(xp.ones(coords0.shape[0], dtype="float32"))
+    tof_backproj = proj.adjoint(xp.ones(coords0.shape[0], dtype="float32"))
 
-xp.save(Path(fname).with_suffix(".tof_backproj.npy"), tof_backproj)
+    tof_backproj_path = Path(fname).with_name(
+        f"{Path(fname).stem}_tof_backproj_{params_hash}.npy"
+    )
+    xp.save(tof_backproj_path, tof_backproj)
 
 del proj
 
@@ -412,7 +454,7 @@
         )
         recon *= tmp / sens_img
 
-opath = Path(fname).parent / f"lm_osem_{num_epochs}_{num_subsets}.npy"
+opath = Path(fname).with_name(f"lm_osem_{num_epochs}_{num_subsets}_{params_hash}.npy")
 xp.save(opath, recon)
 print(f"LM OSEM recon saved to {opath}")