Feature/photons manual calculation

.env_template

@@ -0,0 +1,29 @@
+#                       NOTES
+#
+# 1) A subject-like dataset has the following structure:
+# subject_like_dataset
+# ├── 0
+# .   ├── file0.ext0
+#     .
+#     .
+#     .
+# .   └── file1.ext1
+# 
+# 
+# .
+# └── 99
+#     ├── file0.ext0
+#     .
+#     .
+#     .
+#     └── file1.ext1
+#
+# 2) Change the name of the file to .env or
+# use a symlink
+
+# Env vars for scripts/launch_experiment_set.sh
+export DATA= # path to subject-like activity maps dataset
+export RESULTS= # path to results dir
+export SIMULATIONS= # path to store simulations
+export VENV= # path to python3.9 virtual environment with the requirements.txt installed
+export SIMPET= # path to simpet repository

scripts/experiment.py

@@ -42,16 +42,19 @@ def simulate(cfg: DictConfig) -> None:
 
     OmegaConf.resolve(cfg)
 
-    patient_dirname = cfg.get("params")["patient_dirname"]
     results_dir = Path(cfg["dir_results_path"])
-    patient_dir = results_dir.joinpath(patient_dirname)
+    patient_dirname = str(cfg['params']['patient_dirname'])
+    output_dir = results_dir.joinpath(cfg["params"]["output_dir"])
+    scanner_name = cfg['params']['scanner']["scanner_name"].lower().replace(" ", "_")
+    config_name = f"{patient_dirname}_{scanner_name}.yaml"
+    config_path = output_dir.joinpath(config_name)
 
     test = simpet.SimPET(cfg)
     test.run()
 
     cfg = OmegaConf.to_container(cfg)
 
-    save_cfg(cfg, patient_dir)
+    save_cfg(cfg, config_path)
 
 
 if __name__ == "__main__":

scripts/launch_experiment_set.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+# Variable are defined in the ENVFILE
+
+SCRIPT_DIR="$(dirname "$(readlink -f "$0")")"
+SIMPET_DIR="$(dirname "$SCRIPT_DIR")"
+LOGS_DIR="$SIMPET_DIR"/logs
+ENVFILE="$SIMPET_DIR"/.env
+
+mkdir -p "$LOGS_DIR"
+source "$ENVFILE"
+nohup $VENV -u "$SIMPET"/scripts/experiment_set.py >"$SIMPET"/logs/simulation_set.log 2>&1 &

src/simset/simset_sim.py

@@ -1,16 +1,51 @@
 # -*- coding: utf-8 -*-
-from os.path import join, dirname, abspath, isdir, basename, exists
-import os, sys
+import os
+import re
 import shutil
-import datetime, time
-from multiprocessing import Process
-from utils import resources as rsc
-from utils import tools
-import yaml
+import time
 import numpy as np
-import src.simset.simset_tools as simset_tools
-import subprocess
 import nibabel as nib
+import warnings
+from multiprocessing import Process
+from pathlib import Path
+from os import PathLike
+from os.path import join, dirname, abspath, exists
+import src.simset.simset_tools as simset_tools
+from utils import tools
+
+
+def read_ws_from_simset_log(simset_log: PathLike) -> float:
+    """
+    Read W and W^2 from SimSet log file.
+
+    Args:
+        simset_log: path to simset log file.
+
+    Returns:
+        W^2/W as ``float``.
+    """
+    simset_log = Path(simset_log)
+
+    scintific_number_regex = "(\d+\.\d+e\+\d+)"
+    w_regex = re.compile(
+        f"Sum of accepted coincidence weights in this simulation = {scintific_number_regex}")
+    w2_regex = re.compile(
+        f"Sum of accepted coincidence squared weights in this simulation = {scintific_number_regex}")
+
+    weights = {"W": None, "W2": None}
+    with open(simset_log, 'r') as slog:
+        for line in slog:
+            w_match = w_regex.match(line)
+            w2_match = w2_regex.match(line)
+            if w_match:
+                weights["W"] = float(w_match.groups()[0])
+            elif w2_match:
+                weights["W2"] = float(w2_match.groups()[0])
+            if weights["W"] is not None and weights["W2"] is not None:
+                quotient = weights["W2"] / weights["W"]
+                if quotient < 1:
+                    warnings.warn("W2/W is less than 1!")
+                return quotient
 
 
 class SimSET_Simulation(object):
@@ -107,20 +142,40 @@ def run_simset_simulation(self, sim_dir):
         my_phg = self.prepare_simset_files(
             sim_dir, act_table_factor, act, sim_photons, sim_time, 0
         )
+        # Copying phg file for posterior analysis
+        sim_phg = Path(sim_dir).joinpath("phg.rec")
+        shutil.copy(
+            sim_phg,
+            sim_phg.with_name("phg_first_sampling_sim.rec")
+        )
+        my_log = join(sim_dir, "simset_s0_init.log")
+
+        print("Running first sampling simulation...")
+        command = "%s/bin/phg %s > %s" % (self.simset_dir, my_phg, my_log)
+        tools.osrun(command, log_file)
+
+        my_phg = self.prepare_simset_files(
+            sim_dir, act_table_factor, act, sim_photons, sim_time, 1
+        )
+        # Copying phg file for posterior analysis
+        sim_phg = Path(sim_dir).joinpath("phg.rec")
+        shutil.copy(
+            sim_phg,
+            sim_phg.with_name("phg_second_sampling_sim.rec")
+        )
         my_log = join(sim_dir, "simset_s0.log")
 
+        print("Running second sampling simulation...")
         command = "%s/bin/phg %s > %s" % (self.simset_dir, my_phg, my_log)
         tools.osrun(command, log_file)
+        w_quotient = read_ws_from_simset_log(my_log)
 
         rec_weight = join(sim_dir, "rec.weight")
         det_hf = join(sim_dir, "det_hf.hist")
         phg_hf = join(sim_dir, "phg_hf.hist")
 
         if self.s_photons != 0 and self.params.get("add_randoms") != 1:
-            if self.photons != 0:
-                sim_photons = self.photons / self.divisions
-            else:
-                sim_photons = self.photons
+            sim_photons = int(self.s_photons * w_quotient)
 
             # Removes counts for preparing for the next simulation
             os.remove(rec_weight)
@@ -134,7 +189,7 @@ def run_simset_simulation(self, sim_dir):
             )
             my_log = join(sim_dir, "simset_s1.log")
 
-            print("Running the sencond simulation with importance sampling...")
+            print("Running full simulation with importance sampling...")
 
             command = "%s/bin/phg %s > %s" % (self.simset_dir, my_phg, my_log)
             tools.osrun(command, log_file)

src/stir/stir_tools.py

@@ -24,6 +24,7 @@ def create_stir_hs_from_detparams(scannerParams, output_file, output_format="Sim
     matrix_size, ring_difference = generate_segments_lists_stir(
         num_rings, num_rings - 1
     )
+    transaxial_crystal_distance = scannerParams.get("transaxial_crystal_size")
 
     scanner_name = scannerParams.get("scanner_name")
 
@@ -117,7 +118,7 @@ def create_stir_hs_from_detparams(scannerParams, output_file, output_format="Sim
         + str(td_bins)
         + "\n"
         + "Energy_resolution := "
-        + str(scannerParams.get("energy_resolution") * 5.11)
+        + str(scannerParams.get("energy_resolution") / 511)
         + "\n"
         + "Reference energy (in keV) := 511\n"
         + "Number of blocks per bucket in transaxial direction := 1\n"
@@ -126,6 +127,23 @@ def create_stir_hs_from_detparams(scannerParams, output_file, output_format="Sim
         + "Number of crystals per block in transaxial direction := 1\n"
         + "Number of crystals per singles unit in axial direction := 1\n"
         + "Number of crystals per singles unit in transaxial direction := 1\n"
+        + "Scanner geometry := Cylindrical"
+        + "\n"
+        + "Distance between crystals in axial direction (cm) := "
+        + str(ring_spacing)
+        + "\n"
+        + "Distance between crystals in transaxial direction (cm) := "
+        + str(transaxial_crystal_distance)
+        + "\n"
+        + "Distance between blocks in transaxial direction (cm) := "
+        + str(transaxial_crystal_distance)
+        + "\n"
+        + "Distance between crystals in axial direction (cm) := "
+        + str(ring_spacing)
+        + "\n"
+        + "Distance between blocks in axial direction (cm) := "
+        + str(ring_spacing)
+        + "\n"
         + "end scanner parameters:=\n"
         + "effective central bin size (cm) := "
         + str(bin_size)