From b4f3d02ea54d060632fcd56609959fd40c5761e7 Mon Sep 17 00:00:00 2001 From: David Date: Mon, 7 Sep 2026 14:59:23 +0200 Subject: [PATCH 1/6] dose rate as a proper CLI ; tests redone --- opengate/bin/dose_rate.py | 232 ++++++++++++++++-- opengate/contrib/dose/doserate.py | 91 ++++++- opengate/sources/base.py | 4 +- .../tests/src/source/test035a_dose_rate.py | 5 +- .../src/source/test035b_dose_rate_vrt.py | 177 ++++++++++--- pyproject.toml | 1 + 6 files changed, 438 insertions(+), 72 deletions(-) diff --git a/opengate/bin/dose_rate.py b/opengate/bin/dose_rate.py index af141e7737..83f44aa458 100755 --- a/opengate/bin/dose_rate.py +++ b/opengate/bin/dose_rate.py @@ -1,46 +1,228 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import click +import copy import json +import pathlib +import sys +import click from box import Box -from opengate.contrib.dose.doserate import create_simulation -from opengate.utility import get_random_folder_name +from opengate.contrib.dose.doserate import create_simulation, merge_vrt_dose_rate CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"]) @click.command(context_settings=CONTEXT_SETTINGS) -@click.argument("json_param", nargs=1) -@click.option("--output_folder", "-o", default="AUTO", help="output folder, auto=rnd") -def go(json_param, output_folder): - # open the param file +@click.argument( + "json_param", type=click.Path(exists=True), required=False, default=None +) +@click.option( + "--mode", + "-m", + type=click.Choice( + ["vrt", "analog", "e-", "gamma_tle", "gamma"], case_sensitive=False + ), + default=None, + help="Simulation mode: vrt (runs e- and gamma_tle then merges), analog (full ion decay), e-, gamma_tle, gamma (default: from JSON or 'vrt')", +) +@click.option( + "--activity", + "-a", + default=None, + type=float, + help="Total simulated activity in Bq (overrides JSON). In vrt mode, this applies to gamma.", +) +@click.option( + "--e-factor", + default=10.0, + type=float, + show_default=True, + help="In vrt mode, factor by which electron activity is reduced (activity_e = activity / e_factor).", +) +@click.option( + "--threads", + "-t", + default=None, + type=int, + help="Number of threads (default: from JSON or 4, 1 on Windows).", +) +@click.option( + "--output_folder", + "-o", + default=None, + type=click.Path(), + help="Output folder. Default is auto-named based on mode (e.g. output_vrt, output_analog, output_e, output_gamma_tle).", +) +@click.option( + "--visu", + is_flag=True, + default=False, + help="Enable visualization (forces single thread).", +) +@click.option( + "--merge-only", + nargs=2, + type=click.Path(exists=True), + default=None, + metavar="DIR_E DIR_GAMMA", + help="Merge precomputed electron and gamma simulation folders and exit.", +) +def go(json_param, mode, activity, e_factor, threads, output_folder, visu, merge_only): + # Handle merge-only mode first + if merge_only: + dir_e, dir_gamma = merge_only + out_dir = ( + pathlib.Path(output_folder) + if output_folder + else pathlib.Path("output_merged") + ) + print( + f"Merging VRT dose rate outputs from {dir_e} and {dir_gamma} into {out_dir} (e_factor={e_factor})..." + ) + merged = merge_vrt_dose_rate(dir_e, dir_gamma, out_dir, e_factor=e_factor) + print(f"Merged {len(merged)} image(s):") + for f in merged: + print(f" {f}") + return + + if not json_param: + click.echo( + "Error: Missing argument 'JSON_PARAM'. Use --help for usage instructions.", + err=True, + ) + sys.exit(1) + + # Open the parameter file + json_path = pathlib.Path(json_param).resolve() try: - f = open(json_param, "r") - param = json.load(f) + with open(json_path, "r") as f: + param_dict = json.load(f) except IOError: - print(f"Cannot open input json file {json_param}") - param = Box(param) - print(param) + click.echo(f"Cannot open input json file {json_param}", err=True) + sys.exit(1) - # set or create output_folder - if output_folder == "AUTO": - output_folder = get_random_folder_name() - param.output_folder = output_folder + param = Box(param_dict) - # set activity as int (to deal with 1e4 notation) + # Resolve relative paths in JSON relative to json file directory if not found in cwd + json_dir = json_path.parent + for key in ["ct_image", "table_mat", "table_density", "activity_image"]: + if hasattr(param, key) and param[key]: + p = pathlib.Path(param[key]) + if not p.exists() and (json_dir / p).exists(): + param[key] = str((json_dir / p).resolve()) + + # Apply overrides from CLI + if activity is not None: + param.activity_bq = activity + elif not hasattr(param, "activity_bq") or param.activity_bq is None: + param.activity_bq = 1e6 param.activity_bq = int(float(param.activity_bq)) - # create the simu - sim = create_simulation(param) + default_threads = 1 if sys.platform.startswith("win") else 4 + if threads is not None: + param.number_of_threads = threads + elif not hasattr(param, "number_of_threads") or param.number_of_threads is None: + param.number_of_threads = default_threads + + if visu: + param.visu = True + elif not hasattr(param, "visu"): + param.visu = False + + if not hasattr(param, "verbose"): + param.verbose = True + if not hasattr(param, "density_tolerance_gcm3"): + param.density_tolerance_gcm3 = 0.05 + + # Determine mode: CLI takes precedence, then JSON, then default to "vrt" + if mode is None: + if hasattr(param, "mode") and param.mode: + mode = param.mode + else: + mode = "vrt" + if mode == "": + mode = "analog" + + # Determine output folder + if output_folder: + out_dir = pathlib.Path(output_folder) + elif ( + hasattr(param, "output_folder") + and param.output_folder + and param.output_folder != "AUTO" + ): + out_dir = pathlib.Path(param.output_folder) + else: + out_dir = pathlib.Path(f"output_{mode}") + + # Mode: VRT (automated pipeline) + if mode == "vrt": + print(f"=== Running VRT Dose Rate Simulation ===") + print(f"Target activity (gamma): {param.activity_bq:.2e} Bq") + act_e = int(float(param.activity_bq) / e_factor) + print(f"Electron activity (e_factor={e_factor}): {act_e:.2e} Bq") + print(f"Output directory: {out_dir}") + + out_e = out_dir / "e" + out_gamma = out_dir / "gamma_tle" + out_e.mkdir(parents=True, exist_ok=True) + out_gamma.mkdir(parents=True, exist_ok=True) + + # 1. Run electrons + print("\n--- Step 1/2: Simulating electrons (e-) ---") + param_e = copy.deepcopy(param) + param_e.mode = "e-" + param_e.activity_bq = act_e + param_e.output_folder = str(out_e) + sim_e = create_simulation(param_e) + sim_e.run(start_new_process=True) + stats_e = sim_e.get_actor("Stats") + if stats_e: + print(stats_e) + + # 2. Run photons with TLE + print("\n--- Step 2/2: Simulating photons (gamma_tle) ---") + param_gamma = copy.deepcopy(param) + param_gamma.mode = "gamma_tle" + param_gamma.activity_bq = param.activity_bq + param_gamma.output_folder = str(out_gamma) + sim_gamma = create_simulation(param_gamma) + sim_gamma.run(start_new_process=True) + stats_gamma = sim_gamma.get_actor("Stats") + if stats_gamma: + print(stats_gamma) + + # 3. Merge + print(f"\n--- Merging VRT outputs into {out_dir} (e_factor={e_factor}) ---") + merged = merge_vrt_dose_rate(out_e, out_gamma, out_dir, e_factor=e_factor) + print(f"Merged {len(merged)} image(s):") + for f in merged: + print(f" {f}") + print(f"\nDone. VRT results saved in: {out_dir}") + + else: + # Single simulation mode: analog, e-, gamma_tle, gamma + mode_map = { + "analog": "", + "e-": "e-", + "gamma_tle": "gamma_tle", + "gamma": "gamma", + } + param.mode = mode_map[mode] + out_dir.mkdir(parents=True, exist_ok=True) + param.output_folder = str(out_dir) + + print(f"=== Running Dose Rate Simulation (mode={mode}) ===") + print(f"Simulated activity: {param.activity_bq:.2e} Bq") + print(f"Output directory: {out_dir}") - # run - sim.run() + sim = create_simulation(param) + sim.run(start_new_process=True) - # print results at the end - stats = sim.get_actor("Stats") - print(stats) - print(f"Output in {param.output_folder}") + stats = sim.get_actor("Stats") + if stats: + print(stats) + print(f"\nDone. Results saved in: {out_dir}") # -------------------------------------------------------------------------- diff --git a/opengate/contrib/dose/doserate.py b/opengate/contrib/dose/doserate.py index b33dfa5c98..8d919ea9d3 100755 --- a/opengate/contrib/dose/doserate.py +++ b/opengate/contrib/dose/doserate.py @@ -2,7 +2,6 @@ # -*- coding: utf-8 -*- import pathlib - from opengate.geometry.materials import HounsfieldUnit_to_material from opengate.image import get_translation_between_images_center, read_image_info from opengate.logger import INFO @@ -88,7 +87,7 @@ def create_simulation(param): source.particle = "ion" source.ion.Z = rad_list[param.radionuclide]["Z"] source.ion.A = rad_list[param.radionuclide]["A"] - source.activity = param.activity_bq * Bq / sim.number_of_threads + source.activity = param.activity_bq * Bq source.image = param.activity_image source.direction.type = "iso" source.energy.mono = 0 * keV @@ -127,7 +126,7 @@ def create_simulation(param): print(f"tle_threshold = {dose.tle_threshold/keV} keV") else: dose = sim.add_actor("DoseActor", "dose") - dose.output_filename = "edep.mhd" + dose.output_filename = "output.mhd" dose.dose_uncertainty.active = True dose.dose_squared.active = True dose.dose.active = True @@ -151,3 +150,89 @@ def create_simulation(param): stats.stats.write_to_disk = True return sim + + +def merge_vrt_dose_rate(folder_e, folder_gamma, output_folder, e_factor=1.0): + """ + Merge the dose rate output images from electron and gamma simulations. + Dose_total = e_factor * Dose_e + Dose_gamma + """ + import itk + import numpy as np + import shutil + + folder_e = pathlib.Path(folder_e) + folder_gamma = pathlib.Path(folder_gamma) + output_folder = pathlib.Path(output_folder) + output_folder.mkdir(parents=True, exist_ok=True) + + merged_files = [] + # Determine filename prefix: prefer 'output', fallback to 'edep' + prefix = "output" + if ( + not (folder_e / "output_edep.mhd").exists() + and (folder_e / "edep_edep.mhd").exists() + ): + prefix = "edep" + + # Merge edep and dose images + for suffix in ["edep.mhd", "dose.mhd"]: + image_name = f"{prefix}_{suffix}" + path_e = folder_e / image_name + path_gamma = folder_gamma / image_name + if path_e.exists() and path_gamma.exists(): + img_e = itk.imread(str(path_e)) + img_gamma = itk.imread(str(path_gamma)) + arr_e = itk.GetArrayFromImage(img_e) + arr_gamma = itk.GetArrayFromImage(img_gamma) + + arr_total = float(e_factor) * arr_e + arr_gamma + img_total = itk.GetImageFromArray(arr_total) + img_total.CopyInformation(img_e) + + out_path = output_folder / f"output_{suffix}" + itk.imwrite(img_total, str(out_path)) + merged_files.append(out_path) + + # Merge relative uncertainty if available + unc_name = f"{prefix}_dose_uncertainty.mhd" + dose_name = f"{prefix}_dose.mhd" + unc_e_path = folder_e / unc_name + unc_gamma_path = folder_gamma / unc_name + dose_e_path = folder_e / dose_name + dose_gamma_path = folder_gamma / dose_name + + if ( + unc_e_path.exists() + and unc_gamma_path.exists() + and dose_e_path.exists() + and dose_gamma_path.exists() + ): + arr_unc_e = itk.GetArrayFromImage(itk.imread(str(unc_e_path))) + arr_unc_gamma = itk.GetArrayFromImage(itk.imread(str(unc_gamma_path))) + arr_dose_e = itk.GetArrayFromImage(itk.imread(str(dose_e_path))) + arr_dose_gamma = itk.GetArrayFromImage(itk.imread(str(dose_gamma_path))) + + arr_dose_total = float(e_factor) * arr_dose_e + arr_dose_gamma + # Absolute variance: (F * unc_e * dose_e)^2 + (unc_gamma * dose_gamma)^2 + abs_var = (float(e_factor) * arr_unc_e * arr_dose_e) ** 2 + ( + arr_unc_gamma * arr_dose_gamma + ) ** 2 + with np.errstate(divide="ignore", invalid="ignore"): + arr_unc_total = np.where( + arr_dose_total > 0, np.sqrt(abs_var) / arr_dose_total, 0.0 + ) + + img_unc_total = itk.GetImageFromArray(arr_unc_total.astype(np.float32)) + img_unc_total.CopyInformation(itk.imread(str(unc_e_path))) + out_unc_path = output_folder / "output_dose_uncertainty.mhd" + itk.imwrite(img_unc_total, str(out_unc_path)) + merged_files.append(out_unc_path) + + # Copy labels from electron run if present + for label_file in ["labels.mhd", "labels.raw", "labels.json"]: + src = folder_e / label_file + if src.exists(): + shutil.copy2(src, output_folder / label_file) + + return merged_files diff --git a/opengate/sources/base.py b/opengate/sources/base.py index bed8fe450e..e14dc59b96 100644 --- a/opengate/sources/base.py +++ b/opengate/sources/base.py @@ -246,8 +246,8 @@ def gather_outputs(self, thread_sources): pass def recover_user_output(self, s): - pid = os.getpid() - print(f"(python) recover_user_output {self.name} pid={pid}") + # pid = os.getpid() + # print(f"(python) recover_user_output {self.name} pid={pid}") for k, v in s.user_info.items(): self.user_info[k] = v diff --git a/opengate/tests/src/source/test035a_dose_rate.py b/opengate/tests/src/source/test035a_dose_rate.py index 896e369537..b7eac05850 100755 --- a/opengate/tests/src/source/test035a_dose_rate.py +++ b/opengate/tests/src/source/test035a_dose_rate.py @@ -44,9 +44,6 @@ MeV = gate.g4_units.MeV source.energy.mono = 1 * MeV - print("Phys list cuts:") - print(sim.physics_manager.dump_production_cuts()) - # run sim.run(start_new_process=True) @@ -63,7 +60,7 @@ print(h) is_ok = ( utility.assert_images( - paths.output_ref / "edep.mhd", + paths.output_ref / " .mhd", h.edep.get_output_path(), stats, tolerance=15, diff --git a/opengate/tests/src/source/test035b_dose_rate_vrt.py b/opengate/tests/src/source/test035b_dose_rate_vrt.py index 193e77116d..9c6f3ef443 100755 --- a/opengate/tests/src/source/test035b_dose_rate_vrt.py +++ b/opengate/tests/src/source/test035b_dose_rate_vrt.py @@ -2,76 +2,177 @@ # -*- coding: utf-8 -*- from pathlib import Path - +import json +import sys +import numpy as np import itk from box import Box - import opengate as gate -from opengate.contrib.dose.doserate import create_simulation +from opengate.contrib.dose.doserate import create_simulation, merge_vrt_dose_rate from opengate.tests import utility -def run_simulation_vrt(vrt_mode=""): - paths = utility.get_default_test_paths( - __file__, "", output_folder="test035b_" + vrt_mode - ) +def run_simulation(vrt_mode="analog", activity=5e5, threads=None, skip=False): + if threads is None: + threads = 1 if sys.platform.startswith("win") else 4 + folder_name = "test035b_" + vrt_mode + paths = utility.get_default_test_paths(__file__, "", output_folder=folder_name) dr_data = paths.data / "dose_rate_data" - # set param - gcm3 = gate.g4_units.g_cm3 + # Set parameters param = Box() param.ct_image = str(dr_data / "29_CT_5mm_crop.mhd") param.table_mat = str(dr_data / "Schneider2000MaterialsTable.txt") param.table_density = str(dr_data / "Schneider2000DensitiesTable.txt") param.activity_image = str(dr_data / "activity_test_crop_4mm.mhd") param.radionuclide = "Lu177" - param.activity_bq = 5e5 - param.number_of_threads = 1 + param.activity_bq = activity + param.number_of_threads = threads param.visu = False param.verbose = True param.density_tolerance_gcm3 = 0.05 param.output_folder = str(paths.output) - param.mode = vrt_mode + param.mode = "" if vrt_mode == "analog" else vrt_mode - # Create the simu - # Note that the returned sim object can be modified to change source or cuts or whatever other parameters + # Create the simulation sim = create_simulation(param) - # stats + # Stats actor stats = sim.get_actor("Stats") - stats.output_filename = "stats035_" + vrt_mode + ".txt" - - print("Phys list cuts:") - print(sim.physics_manager.dump_production_cuts()) + stats.output_filename = f"stats035b_{vrt_mode}.txt" - # run - if not vrt_mode == "": + # Run in a new process + if not skip: sim.run(start_new_process=True) + return paths if __name__ == "__main__": - paths_original_simu = run_simulation_vrt( - vrt_mode="" - ) # too long do not run the simulation without vrt, use reference output - paths_vrt_e_simu = run_simulation_vrt(vrt_mode="e-") - paths_vrt_gamma_simu = run_simulation_vrt(vrt_mode="gamma") + activity = 5e5 + e_factor = 1 + + # 1. Analog simulation + print("Analog simulation ...") + paths_analog = run_simulation(vrt_mode="analog", activity=activity, skip=False) + analog_edep = paths_analog.output / "output_edep.mhd" + analog_unc = paths_analog.output / "output_dose_uncertainty.mhd" + + # 2. VRT simulations using the new API (e- and gamma_tle) + print() + print("VRT simulations (e-, low stats) ..") + paths_vrt_e = run_simulation( + vrt_mode="e-", activity=int(activity / e_factor), skip=False + ) - # dose comparison between original simu and vrt simu (electron + gamma) print() - gate.exception.warning(f"Check dose") - dose_vrt_e_simu = itk.imread(paths_vrt_e_simu.output / "edep_edep.mhd") - dose_vrt_gamma_simu = itk.imread(paths_vrt_gamma_simu.output / "edep_edep.mhd") - array_vrt_e_simu = itk.GetArrayFromImage(dose_vrt_e_simu) - array_vrt_gamma_simu = itk.GetArrayFromImage(dose_vrt_gamma_simu) - array_vrt = array_vrt_e_simu + array_vrt_gamma_simu - dose_vrt = itk.GetImageFromArray(array_vrt) - dose_vrt.CopyInformation(dose_vrt_e_simu) - itk.imwrite(dose_vrt, paths_original_simu.output / "edep_edep_vrt.mhd") + print("VRT simulations (gamma tle) ..") + paths_vrt_gamma = run_simulation(vrt_mode="gamma_tle", activity=activity) + + # 3. Merge outputs using merge_vrt_dose_rate + paths_vrt_merged = utility.get_default_test_paths( + __file__, "", output_folder="test035b_vrt" + ) + print( + f"\nMerging VRT outputs into {paths_vrt_merged.output} (e_factor={e_factor})..." + ) + merged_files = merge_vrt_dose_rate( + paths_vrt_e.output, + paths_vrt_gamma.output, + paths_vrt_merged.output, + e_factor=e_factor, + ) + print(f"Merged {len(merged_files)} files: {[str(f) for f in merged_files]}") + + vrt_edep = paths_vrt_merged.output / "output_edep.mhd" + vrt_dose = paths_vrt_merged.output / "output_dose.mhd" + vrt_unc = paths_vrt_merged.output / "output_dose_uncertainty.mhd" + + # 4. Check dose (edep) + print() + gate.exception.warning("Check dose (edep)") is_ok = utility.assert_images( - paths_original_simu.output_ref / "edep_edep.mhd", - paths_original_simu.output / "edep_edep_vrt.mhd", + analog_edep, + vrt_edep, tolerance=30, ignore_value_data2=0, ) + + # 5. Check uncertainty & efficiency + print() + gate.exception.warning("Check dose uncertainty & efficiency") + if vrt_unc.exists(): + img_unc_vrt = itk.imread(str(vrt_unc)) + arr_unc_vrt = itk.GetArrayFromImage(img_unc_vrt) + arr_dose_vrt = itk.GetArrayFromImage(itk.imread(str(vrt_dose))) + mask_vrt = arr_dose_vrt > 0 + + mean_unc_vrt = ( + float(np.mean(arr_unc_vrt[mask_vrt])) if np.any(mask_vrt) else 0.0 + ) + var_vrt = ( + float(np.mean(arr_unc_vrt[mask_vrt] ** 2)) if np.any(mask_vrt) else 0.0 + ) + unc_valid = (0.0 < mean_unc_vrt < 1.0) and not np.isnan(mean_unc_vrt) + utility.print_test( + unc_valid, f"VRT dose uncertainty is valid (mean = {mean_unc_vrt:.4f})" + ) + is_ok = is_ok and unc_valid + + if analog_unc.exists(): + img_unc_analog = itk.imread(str(analog_unc)) + arr_unc_analog = itk.GetArrayFromImage(img_unc_analog) + arr_dose_analog = itk.GetArrayFromImage(itk.imread(str(analog_edep))) + mask_analog = arr_dose_analog > 0 + + mean_unc_analog = ( + float(np.mean(arr_unc_analog[mask_analog])) + if np.any(mask_analog) + else 0.0 + ) + var_analog = ( + float(np.mean(arr_unc_analog[mask_analog] ** 2)) + if np.any(mask_analog) + else 0.0 + ) + + print(f"Analog dose uncertainty: mean = {mean_unc_analog:.4f}") + print(f"VRT dose uncertainty: mean = {mean_unc_vrt:.4f}") + + # Variance reduction check: VRT uncertainty should be lower than analog + is_reduced = mean_unc_vrt < mean_unc_analog + utility.print_test( + is_reduced, + f"VRT uncertainty is lower than analog ({mean_unc_vrt:.4f} < {mean_unc_analog:.4f})", + ) + is_ok = is_ok and is_reduced + + # Compute efficiency & speedup + try: + stats_a = paths_analog.output / "stats035b_analog.txt" + stats_e = paths_vrt_e.output / "stats035b_e-.txt" + stats_g = paths_vrt_gamma.output / "stats035b_gamma_tle.txt" + if stats_a.exists() and stats_e.exists() and stats_g.exists(): + with open(stats_a) as f: + t_analog = json.load(f)["duration"]["value"] + with open(stats_e) as f: + t_e = json.load(f)["duration"]["value"] + with open(stats_g) as f: + t_g = json.load(f)["duration"]["value"] + t_vrt = t_e + t_g + + eff_analog = 1.0 / (t_analog**2 * var_analog) + eff_vrt = 1.0 / (t_vrt**2 * var_vrt) + speedup = eff_vrt / eff_analog + + print( + f"\nTime: Analog = {t_analog:.2f} s, VRT = {t_vrt:.2f} s (e-: {t_e:.2f} s, gamma: {t_g:.2f} s)" + ) + print( + f"Efficiency 1 / (t^2 * var): Analog = {eff_analog:.4e}, VRT = {eff_vrt:.4e}" + ) + print(f"Estimated speedup: {speedup:.2f}x") + except Exception as e: + print(f"Could not compute speedup: {e}") + utility.test_ok(is_ok) diff --git a/pyproject.toml b/pyproject.toml index 139dc7912a..fed9a0bd32 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,6 +56,7 @@ opengate_jobs_status = "opengate.bin.opengate_jobs_status:go" opengate_jobs_merge = "opengate.bin.opengate_jobs_merge:go" opengate_jobs_clean = "opengate.bin.opengate_jobs_clean:go" opengate_job_runner = "opengate.bin.opengate_job_runner:go" +opengate_dose_rate = "opengate.bin.dose_rate:go" dose_rate = "opengate.bin.dose_rate:go" split_spect_projections = "opengate.bin.split_spect_projections:go" From b09e45d064173bd309112bca871a0cf578e96db5 Mon Sep 17 00:00:00 2001 From: David Date: Wed, 9 Sep 2026 15:48:16 +0200 Subject: [PATCH 2/6] improve dose rate cli --- opengate/bin/dose_rate.py | 69 ++++++++++++++++++++--- opengate/contrib/dose/doserate.py | 91 +++++++++++++++++++++++++++---- 2 files changed, 141 insertions(+), 19 deletions(-) diff --git a/opengate/bin/dose_rate.py b/opengate/bin/dose_rate.py index 83f44aa458..bc4cbd748b 100755 --- a/opengate/bin/dose_rate.py +++ b/opengate/bin/dose_rate.py @@ -32,12 +32,18 @@ type=float, help="Total simulated activity in Bq (overrides JSON). In vrt mode, this applies to gamma.", ) +@click.option( + "--radionuclide", + "-r", + default=None, + type=str, + help="Radionuclide name (e.g. Lu177, Y90, Ac225) or ion 'Z A' (e.g. '89 225') (overrides JSON).", +) @click.option( "--e-factor", - default=10.0, + default=None, type=float, - show_default=True, - help="In vrt mode, factor by which electron activity is reduced (activity_e = activity / e_factor).", + help="In vrt mode, factor by which electron activity is reduced (default: from JSON or 10.0).", ) @click.option( "--threads", @@ -67,7 +73,17 @@ metavar="DIR_E DIR_GAMMA", help="Merge precomputed electron and gamma simulation folders and exit.", ) -def go(json_param, mode, activity, e_factor, threads, output_folder, visu, merge_only): +def go( + json_param, + mode, + activity, + radionuclide, + e_factor, + threads, + output_folder, + visu, + merge_only, +): # Handle merge-only mode first if merge_only: dir_e, dir_gamma = merge_only @@ -76,10 +92,11 @@ def go(json_param, mode, activity, e_factor, threads, output_folder, visu, merge if output_folder else pathlib.Path("output_merged") ) + ef = e_factor if e_factor is not None else 10.0 print( - f"Merging VRT dose rate outputs from {dir_e} and {dir_gamma} into {out_dir} (e_factor={e_factor})..." + f"Merging VRT dose rate outputs from {dir_e} and {dir_gamma} into {out_dir} (e_factor={ef})..." ) - merged = merge_vrt_dose_rate(dir_e, dir_gamma, out_dir, e_factor=e_factor) + merged = merge_vrt_dose_rate(dir_e, dir_gamma, out_dir, e_factor=ef) print(f"Merged {len(merged)} image(s):") for f in merged: print(f" {f}") @@ -112,6 +129,9 @@ def go(json_param, mode, activity, e_factor, threads, output_folder, visu, merge param[key] = str((json_dir / p).resolve()) # Apply overrides from CLI + if radionuclide is not None: + param.radionuclide = radionuclide + if activity is not None: param.activity_bq = activity elif not hasattr(param, "activity_bq") or param.activity_bq is None: @@ -129,20 +149,55 @@ def go(json_param, mode, activity, e_factor, threads, output_folder, visu, merge elif not hasattr(param, "visu"): param.visu = False + if e_factor is not None: + param.e_factor = e_factor + elif "e_factor" in param and param.e_factor is not None: + param.e_factor = float(param.e_factor) + elif "e-factor" in param and param["e-factor"] is not None: + param.e_factor = float(param["e-factor"]) + else: + param.e_factor = 10.0 + e_factor = float(param.e_factor) + if not hasattr(param, "verbose"): param.verbose = True if not hasattr(param, "density_tolerance_gcm3"): param.density_tolerance_gcm3 = 0.05 - # Determine mode: CLI takes precedence, then JSON, then default to "vrt" + # Determine if radionuclide is an alpha emitter or generic ion 'Z A' (incompatible with VRT) + is_non_vrt_rad = False + if hasattr(param, "radionuclide"): + rad_val = param.radionuclide + if isinstance(rad_val, (list, tuple)): + is_non_vrt_rad = True + elif isinstance(rad_val, str): + rad_str = rad_val.strip() + if rad_str in ["Ac225", "Ra223", "Bi213", "Pb212"]: + is_non_vrt_rad = True + elif rad_str.startswith("ion") or ( + len(rad_str.split()) > 1 and rad_str.split()[0].isdigit() + ): + is_non_vrt_rad = True + + # Determine mode: CLI takes precedence, then JSON, then default if mode is None: if hasattr(param, "mode") and param.mode: mode = param.mode + elif is_non_vrt_rad: + mode = "analog" else: mode = "vrt" if mode == "": mode = "analog" + if mode == "vrt" and is_non_vrt_rad: + click.echo( + f"Error: VRT mode is designed for beta/gamma emitters (e.g. Lu177, Y90). " + f"For '{param.radionuclide}', please use --mode analog.", + err=True, + ) + sys.exit(1) + # Determine output folder if output_folder: out_dir = pathlib.Path(output_folder) diff --git a/opengate/contrib/dose/doserate.py b/opengate/contrib/dose/doserate.py index 8d919ea9d3..a117a931f2 100755 --- a/opengate/contrib/dose/doserate.py +++ b/opengate/contrib/dose/doserate.py @@ -10,6 +10,80 @@ from opengate.utility import g4_best_unit, g4_units +def get_ion_z_a_e(rad): + """ + Parse a radionuclide specification into (Z, A, E). + Accepts: + - Dictionary lookup: 'Lu177', 'Y90', 'In111', 'I131', 'Ac225', etc. + - Numeric format: '89 225', 'ion 89 225', '89, 225', '89-225', [89, 225] + - Element notation: 'Ac225', 'Ac-225', '225Ac', 'Tb161' + """ + rad_list = { + "Lu177": {"Z": 71, "A": 177}, + "Y90": {"Z": 39, "A": 90}, + "In111": {"Z": 49, "A": 111}, + "I131": {"Z": 53, "A": 131}, + "Ac225": {"Z": 89, "A": 225}, + "Ra223": {"Z": 88, "A": 223}, + "Bi213": {"Z": 83, "A": 213}, + "Pb212": {"Z": 82, "A": 212}, + "Tb161": {"Z": 65, "A": 161}, + } + if isinstance(rad, (list, tuple)): + z = int(rad[0]) + a = int(rad[1]) + e = int(rad[2]) if len(rad) > 2 else 0 + return z, a, e + + if not isinstance(rad, str): + raise ValueError(f"Unsupported radionuclide specification: {rad}") + + rad_str = rad.strip() + if rad_str in rad_list: + return rad_list[rad_str]["Z"], rad_list[rad_str]["A"], 0 + + if rad_str.startswith("ion"): + parts = rad_str.split() + z = int(parts[1]) + a = int(parts[2]) + e = int(parts[3]) if len(parts) > 3 else 0 + return z, a, e + + # Format: "89 225" or "89, 225" or "89-225" + import re + + m = re.match(r"^(\d+)[,\s_-]+(\d+)(?:[,\s_-]+(\d+))?$", rad_str) + if m: + z = int(m.group(1)) + a = int(m.group(2)) + e = int(m.group(3)) if m.group(3) else 0 + return z, a, e + + import opengate_core as g4 + + # Element symbol + mass: "Ac225", "Ac-225" + m = re.match(r"^([a-zA-Z]+)[-_]?(\d+)$", rad_str) + if m: + elem = m.group(1).capitalize() + a = int(m.group(2)) + z = g4.G4NistManager.Instance().GetZ(elem) + if z == 0: + raise ValueError(f"Unknown element symbol '{elem}' in radionuclide '{rad}'") + return z, a, 0 + + # Mass + element symbol: "225Ac" + m = re.match(r"^(\d+)[-_]?([a-zA-Z]+)$", rad_str) + if m: + a = int(m.group(1)) + elem = m.group(2).capitalize() + z = g4.G4NistManager.Instance().GetZ(elem) + if z == 0: + raise ValueError(f"Unknown element symbol '{elem}' in radionuclide '{rad}'") + return z, a, 0 + + raise ValueError(f"Cannot parse radionuclide name or ion definition '{rad}'") + + def create_simulation(param): """ param is dict with: @@ -71,22 +145,15 @@ def create_simulation(param): ) ct.dump_label_image = param.output_folder / "labels.mhd" - # some radionuclides choice - # (user of this function can still change - # the source in the output sim) - rad_list = { - "Lu177": {"Z": 71, "A": 177, "name": "Lutetium 177"}, - "Y90": {"Z": 39, "A": 90, "name": "Yttrium 90"}, - "In111": {"Z": 49, "A": 111, "name": "Indium 111"}, - "I131": {"Z": 53, "A": 131, "name": "Iodine 131"}, - } - # Activity source from an image source = sim.add_source("VoxelSource", "vox") source.attached_to = ct.name source.particle = "ion" - source.ion.Z = rad_list[param.radionuclide]["Z"] - source.ion.A = rad_list[param.radionuclide]["A"] + z, a, e = get_ion_z_a_e(param.radionuclide) + source.ion.Z = z + source.ion.A = a + if e: + source.ion.E = e source.activity = param.activity_bq * Bq source.image = param.activity_image source.direction.type = "iso" From b12e225c5ae166bd6dcf61dc692a053a27f028fc Mon Sep 17 00:00:00 2001 From: David Date: Wed, 9 Sep 2026 15:56:17 +0200 Subject: [PATCH 3/6] update doc --- .../user_guide_contrib_dose_rate.rst | 253 +++++++++++------- 1 file changed, 158 insertions(+), 95 deletions(-) diff --git a/docs/source/user_guide/user_guide_contrib_dose_rate.rst b/docs/source/user_guide/user_guide_contrib_dose_rate.rst index 38226b24ee..10756a9c70 100644 --- a/docs/source/user_guide/user_guide_contrib_dose_rate.rst +++ b/docs/source/user_guide/user_guide_contrib_dose_rate.rst @@ -4,83 +4,160 @@ Dose rate computation ===================== Dose rate computations can be performed using Monte Carlo simulations, especially in the context of internal dosimetry for targeted radionuclide therapy (TRT). +OpenGATE provides dedicated tools in `opengate.contrib.dose.doserate `_ and the command-line interface ``opengate_dose_rate``. Command-line usage ------------------ -To run the simulation from the command line, use the ``dose_rate`` executable with a JSON configuration file: +To run the simulation from the command line, use the ``opengate_dose_rate`` executable (also available as ``dose_rate``) with a JSON configuration file: .. code-block:: bash - dose_rate dose_rate_test1.json -o outputFolder/ - -The ``dose_rate_test1.json`` file contains the input parameters for the simulation: CT image, material/density calibration tables, activity image (e.g. from SPECT/PET), radionuclide, simulated activity, number of threads, and visualization settings. + opengate_dose_rate dose_rate_param.json + +CLI options +~~~~~~~~~~~ + +.. code-block:: text + + Usage: opengate_dose_rate [OPTIONS] [JSON_PARAM] + + Options: + -m, --mode [vrt|analog|e-|gamma_tle|gamma] + Simulation mode: vrt (runs e- and gamma_tle + then merges), analog (full ion decay), e-, + gamma_tle, gamma (default: from JSON or + 'vrt') + -a, --activity FLOAT Total simulated activity in Bq (overrides + JSON). In vrt mode, this applies to gamma. + -r, --radionuclide TEXT Radionuclide name (e.g. Lu177, Y90, Ac225) + or ion 'Z A' (e.g. '89 225') (overrides JSON). + --e-factor FLOAT In vrt mode, factor by which electron + activity is reduced (default: from JSON or + 10.0). + -t, --threads INTEGER Number of threads (default: from JSON or 4, + 1 on Windows). + -o, --output_folder PATH Output folder. Default is auto-named based + on mode (e.g. output_vrt, output_analog, + output_e, output_gamma_tle). + --visu Enable visualization (forces single thread). + --merge-only DIR_E DIR_GAMMA Merge precomputed electron and gamma + simulation folders and exit. + -h, --help Show this message and exit. + +JSON configuration +~~~~~~~~~~~~~~~~~~ + +The JSON file contains the input parameters for the simulation: .. code-block:: json { - "# Input CT image": "", "ct_image": "./dose_rate_data/29_CT_5mm.mhd", "table_mat": "./dose_rate_data/Schneider2000MaterialsTable.txt", "table_density": "./dose_rate_data/Schneider2000DensitiesTable.txt", - "density_tolerance_gcm3": 0.2, + "density_tolerance_gcm3": 0.05, - "# Input activity image": "", "activity_image": "./dose_rate_data/385_NM_5mm.mhd", - - "# Input radionuclide": "", "radionuclide": "Lu177", - - "# Input total simulated activity in the whole image, in Bq": "", "activity_bq": 1e6, - "# Option: number of threads": "", - "number_of_threads": 1, - - "# Option: visualisation (for debug)": "", + "mode": "vrt", + "e_factor": 10.0, + "number_of_threads": 4, "visu": false, - - "# verbosity": "", "verbose": true } -The SPECT/PET and CT images provide information on the bio-distribution of the source as well as the density and composition of patient tissues. The simulated activity can be lower than the actual administered activity to reduce computation times; the resulting dose rates can subsequently be scaled to the injected activity. -The SPECT/PET and CT images provide information on the bio-distribution of the source as well as the density and composition of patient tissues. Note that ``activity_image`` is used as a **relative spatial probability distribution**: it is internally normalized such that the sum of all voxel values is 1. The absolute total activity (in Bq) simulated across the entire volume is set independently by ``param.activity_bq`` (or ``source.activity``). +Parameter descriptions: -The simulated activity can be lower than the actual administered activity to reduce computation times; the resulting dose rates can subsequently be scaled linearly to the injected activity. +- **ct_image**: Patient CT image providing the geometry and tissue distribution. +- **table_mat**, **table_density**: Calibration tables mapping Hounsfield Units (HU) to materials and mass densities (e.g., Schneider 2000). +- **density_tolerance_gcm3**: Tolerance for grouping voxels with similar densities into discrete material definitions (default: 0.05 g/cm³). +- **activity_image**: 3D spatial activity distribution (from SPECT or PET). Used as a **relative spatial probability distribution** (internally normalized so that the sum of all voxel values is 1). +- **radionuclide**: Radionuclide identifier (see below). +- **activity_bq**: Total simulated activity across the whole volume in Becquerels (Bq). +- **mode**: Simulation mode: + - ``"vrt"``: Runs decoupled electron and gamma-TLE simulations and merges them automatically (default for beta/gamma emitters). + - ``"analog"``: Full ion radioactive decay simulation (default for alpha emitters and custom ions). + - Single-component modes: ``"e-"``, ``"gamma_tle"``, or ``"gamma"``. +- **e_factor** (or ``e-factor``): In VRT mode, the reduction factor for electron histories (``activity_e = activity / e_factor``, default: 10.0). +- **number_of_threads**: Number of worker threads (default: 4 on Linux/macOS, 1 on Windows). +- **visu**: Set to ``true`` to generate a VRML visualization for debugging (forces 1 thread). -By default, the simulation duration is set to 1 second. It can be modified by setting the time intervals corresponding to the acquisition duration: +Radionuclide specification +~~~~~~~~~~~~~~~~~~~~~~~~~~ -.. code-block:: python +The ``radionuclide`` parameter accepts several flexible formats: - sim.run_timing_intervals = [[0, 3600 * sec]] # for a total duration of 1 hour +1. **Standard isotope names**: + - Built-in: ``"Lu177"``, ``"Y90"``, ``"In111"``, ``"I131"``, ``"Ac225"``, ``"Ra223"``, ``"Bi213"``, ``"Pb212"``, ``"Tb161"``. + - Dynamic: Any valid isotope name (e.g. ``"Sm153"``, ``"Re186"``, ``"Cu67"``) is dynamically resolved to its atomic number :math:`Z` and mass :math:`A` using the Geant4 NIST database. -Radioactive decay can also be accounted for by specifying the radionuclide half-life: +2. **Generic ion specification ("Z A")**: + - Space/comma/hyphen separated: ``"89 225"``, ``"ion 89 225"``, ``"89, 225"``, ``"89-225"``. + - List: ``[89, 225]``. + - Optional excitation energy :math:`E` can be provided: ``"ion 89 225 0"``. + This allows simulating any ion without requiring database entries. -.. code-block:: python +.. note:: + - **Beta/gamma emitters** (such as Lu-177, Y-90, I-131, In-111) can run in either accelerated ``vrt`` mode or ``analog`` mode. + - **Alpha emitters** (such as Ac-225, Ra-223) and generic ions (``"89 225"``) do not have a beta spectrum and should be run in ``analog`` mode (``source.particle = "ion"`` with full radioactive decay enabled). When an alpha emitter or generic ion is detected, ``opengate_dose_rate`` automatically selects ``analog`` mode. + +Simulation outputs +~~~~~~~~~~~~~~~~~~ + +The simulation generates the following files in the output directory: + +- **output_dose.mhd**: Absorbed dose rate map (in Gy). +- **output_edep.mhd**: Energy deposition map (in MeV). +- **output_dose_uncertainty.mhd**: Statistical relative uncertainty map (:math:`u = \sigma(D) / D`). +- **labels.mhd**: Material labels segmented from the CT image. +- **stats.txt**: Detailed runtime statistics (histories, tracks, steps, computation duration). + + +Accelerated computation with Variance Reduction Techniques (VRT) +---------------------------------------------------------------- + +In standard analog simulations (``mode = "analog"``), radionuclide decay simulates the parent ion and tracks all decay products (electrons, alphas, gammas, X-rays) in the same process. Because electrons have short ranges and deposit their energy locally, while photons travel long distances with low interaction probability per voxel, tracking both with analog Monte Carlo can be computationally demanding. + +To significantly accelerate dose rate computations for beta/gamma emitters, an automated VRT pipeline is available: + +1. **Emission Decoupling**: + - **Electron simulation** (``mode = "e-"``): Simulates electrons sampled from the radionuclide beta spectrum. A 1 m production cut in the CT volume deposits electron energy locally in the voxel of emission. Because electron dose is deposited locally with high efficiency, far fewer histories are needed (:math:`A_e = A / \text{e\_factor}`, default reduction factor 10). + - **Photon simulation** (``mode = "gamma_tle"``): Simulates photons sampled from the gamma spectrum using a :class:`~.opengate.actors.doseactors.TLEDoseActor` (Track Length Estimator). TLE analytically scores dose along photon paths in every voxel traversed, achieving smooth dose distributions with low uncertainty without tracking rare secondary electrons. - source.half_life = 60 * sec +2. **Automated Merging**: + The ``merge_vrt_dose_rate`` function combines the dose maps: -Supported radionuclides in this helper currently include: ``Lu177``, ``Y90``, ``In111``, and ``I131``. + .. math:: -The simulation generates output files including: -- Dose rate map (``dose_edep.mhd``) -- Energy deposition map (``edep.mhd``) -- Statistical uncertainty map (``edep_uncertainty.mhd``) -- CT material labels (``labels.mhd``) -- Simulation statistics (``stats.txt``) + D_{\text{total}} = F \cdot D_{e} + D_{\gamma} + + where :math:`F` is the ``e_factor``. + Because the electron and gamma simulations are statistically independent, the combined relative uncertainty is computed rigorously via error propagation: + + .. math:: + + u_{\text{total}} = \frac{\sqrt{(F \cdot u_e \cdot D_e)^2 + (u_\gamma \cdot D_\gamma)^2}}{F \cdot D_e + D_\gamma} + + where :math:`u_e` and :math:`u_\gamma` are the relative uncertainties from each run. + +3. **Performance Gain**: + Compared to an analog simulation, the VRT approach with TLE and :math:`F=10` typically achieves an efficiency speedup (:math:`\epsilon = 1 / (t^2 \cdot \text{variance})`) of **10x to 30x** for equivalent statistical precision. Python API usage ---------------- -You can also set up and customize the dose rate simulation directly in Python using the helper module `opengate.contrib.dose.doserate `_ (see ``opengate/tests/src/source/test035a_dose_rate.py``): +Standard simulation (Analog) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +To set up an analog simulation directly in Python (see ``opengate/tests/src/source/test035a_dose_rate.py``): .. code-block:: python - from pathlib import Path from box import Box - import opengate as gate from opengate.contrib.dose.doserate import create_simulation # Configure parameters @@ -89,78 +166,64 @@ You can also set up and customize the dose rate simulation directly in Python us param.table_mat = "Schneider2000MaterialsTable.txt" param.table_density = "Schneider2000DensitiesTable.txt" param.activity_image = "activity_test_crop_4mm.mhd" - param.radionuclide = "Lu177" + param.radionuclide = "Lu177" # or "Ac225", "89 225", etc. param.activity_bq = 1e6 param.number_of_threads = 4 param.visu = False param.verbose = True param.density_tolerance_gcm3 = 0.05 - param.output_folder = "output_test035a" + param.output_folder = "output_analog" param.mode = "" # standard full simulation - # Create the simulation object + # Create and run simulation sim = create_simulation(param) - - # You can customize any component before running sim.run(start_new_process=True) -Accelerated computation with Variance Reduction Techniques (VRT) ----------------------------------------------------------------- - -In standard analog simulations (``param.mode = ""``), radionuclide decay emits both electrons (:math:`\beta^-` / Auger / conversion electrons) and photons (:math:`\gamma` / X-rays) in the same tracking process. Because electrons have short ranges and deposit most of their energy locally, while photons travel long distances with low interaction probability per voxel, tracking both with analog Monte Carlo can be computationally demanding. - -To significantly accelerate dose rate computations, a Variance Reduction Technique (VRT) approach is implemented (see ``opengate/tests/src/source/test035b_dose_rate_vrt.py``): - -1. **Emission Decoupling**: The radionuclide decay is decoupled into separate simulations: - - **Electron simulation** (``param.mode = "e-"``): Emits electrons sampled from the radionuclide beta spectrum. A high production cut (e.g. 1 m in the CT volume) is applied to deposit electron energy locally within the voxel where they originate, drastically speeding up electron scoring. - - **Photon simulation** (``param.mode = "gamma"`` or ``param.mode = "gamma_tle"``): Emits photons sampled from the gamma spectrum. When using ``"gamma_tle"``, the simulation replaces stochastic photon dose deposition with a :class:`~.opengate.actors.doseactors.TLEDoseActor` (Track Length Estimator), achieving faster convergence with fewer simulated particles. +VRT simulation and merging +~~~~~~~~~~~~~~~~~~~~~~~~~~ -2. **Dose Map Recombination**: The dose/energy deposition maps from both runs are summed to obtain the total dose rate. - -Example: running decoupled VRT simulations and combining results -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +To run decoupled VRT simulations and merge their results (see ``opengate/tests/src/source/test035b_dose_rate_vrt.py``): .. code-block:: python - import itk from box import Box - import opengate as gate - from opengate.contrib.dose.doserate import create_simulation - - def run_dose_rate(mode, output_folder, activity_bq=5e5): - param = Box() - param.ct_image = "29_CT_5mm_crop.mhd" - param.table_mat = "Schneider2000MaterialsTable.txt" - param.table_density = "Schneider2000DensitiesTable.txt" - param.activity_image = "activity_test_crop_4mm.mhd" - param.radionuclide = "Lu177" - param.activity_bq = activity_bq - param.number_of_threads = 4 - param.visu = False - param.verbose = True - param.density_tolerance_gcm3 = 0.05 - param.output_folder = output_folder - param.mode = mode # "e-", "gamma", or "gamma_tle" - - sim = create_simulation(param) - sim.run(start_new_process=True) - - # 1. Run electron simulation (local deposition approximation) - run_dose_rate(mode="e-", output_folder="output_vrt_e") - - # 2. Run gamma simulation (TLE or standard photon tracking) - run_dose_rate(mode="gamma", output_folder="output_vrt_gamma") - - # 3. Sum the resulting energy deposition maps - dose_e = itk.imread("output_vrt_e/edep_edep.mhd") - dose_gamma = itk.imread("output_vrt_gamma/edep_edep.mhd") - - array_e = itk.GetArrayFromImage(dose_e) - array_gamma = itk.GetArrayFromImage(dose_gamma) - array_total = array_e + array_gamma - - dose_total = itk.GetImageFromArray(array_total) - dose_total.CopyInformation(dose_e) - itk.imwrite(dose_total, "output_total_vrt/edep_edep_vrt.mhd") - + from opengate.contrib.dose.doserate import create_simulation, merge_vrt_dose_rate + + activity = 1e6 + e_factor = 10.0 + + # Common parameters + base_param = Box() + base_param.ct_image = "29_CT_5mm_crop.mhd" + base_param.table_mat = "Schneider2000MaterialsTable.txt" + base_param.table_density = "Schneider2000DensitiesTable.txt" + base_param.activity_image = "activity_test_crop_4mm.mhd" + base_param.radionuclide = "Lu177" + base_param.number_of_threads = 4 + base_param.density_tolerance_gcm3 = 0.05 + + # 1. Run electron simulation with reduced statistics + param_e = base_param.copy() + param_e.mode = "e-" + param_e.activity_bq = int(activity / e_factor) + param_e.output_folder = "output_e" + sim_e = create_simulation(param_e) + sim_e.run(start_new_process=True) + + # 2. Run photon simulation with TLE + param_g = base_param.copy() + param_g.mode = "gamma_tle" + param_g.activity_bq = activity + param_g.output_folder = "output_gamma_tle" + sim_g = create_simulation(param_g) + sim_g.run(start_new_process=True) + + # 3. Merge outputs and calculate combined uncertainty + merged_files = merge_vrt_dose_rate( + folder_e="output_e", + folder_gamma="output_gamma_tle", + output_folder="output_vrt", + e_factor=e_factor, + ) + print("Merged outputs:", merged_files) From bcbd7ea55cb76656ce48c0b9f6f2e3538f2bb679 Mon Sep 17 00:00:00 2001 From: David Date: Wed, 9 Sep 2026 16:16:45 +0200 Subject: [PATCH 4/6] test for CLI --- .../src/source/test035c_dose_rate_cli.py | 131 ++++++++++++++++++ 1 file changed, 131 insertions(+) create mode 100755 opengate/tests/src/source/test035c_dose_rate_cli.py diff --git a/opengate/tests/src/source/test035c_dose_rate_cli.py b/opengate/tests/src/source/test035c_dose_rate_cli.py new file mode 100755 index 0000000000..70e0e5b0c8 --- /dev/null +++ b/opengate/tests/src/source/test035c_dose_rate_cli.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import json +import pathlib +import tempfile +from click.testing import CliRunner + +from opengate.bin.dose_rate import go +from opengate.tests import utility + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "", output_folder="test035c") + dr_data = paths.data / "dose_rate_data" + + is_ok = True + runner = CliRunner() + + print("=== Test 1: CLI help and argument checking ===") + res_help = runner.invoke(go, ["--help"]) + print(f" --help exit code: {res_help.exit_code}") + if res_help.exit_code != 0 or "--mode" not in res_help.output: + is_ok = False + + res_missing = runner.invoke(go, []) + print(f" Missing json exit code: {res_missing.exit_code}") + if res_missing.exit_code == 0: + is_ok = False + + res_notfound = runner.invoke(go, ["non_existent_file.json"]) + print(f" Non-existent json exit code: {res_notfound.exit_code}") + if res_notfound.exit_code == 0: + is_ok = False + + print("\n=== Test 2: Incompatible mode and radionuclide validation ===") + with tempfile.TemporaryDirectory() as tmpdir: + tmp = pathlib.Path(tmpdir) + dummy_json = tmp / "dummy.json" + dummy_json.write_text( + json.dumps({"ct_image": "dummy.mhd", "radionuclide": "89 225"}) + ) + + # Radionuclide '89 225' (Ac225) with --mode vrt must fail with informative error + res_vrt_fail = runner.invoke(go, [str(dummy_json), "--mode", "vrt"]) + print(f" VRT with ion '89 225' exit code: {res_vrt_fail.exit_code}") + if ( + res_vrt_fail.exit_code == 0 + or "Error: VRT mode is designed" not in res_vrt_fail.output + ): + is_ok = False + else: + print(" Rejection message verified successfully.") + + print("\n=== Test 3: --merge-only mode ===") + dir_e = paths.output_ref.parent / "test035b" / "test035b_e-" + if not dir_e.exists(): + dir_e = paths.output.parent / "test035b_e-" + dir_gamma = paths.output_ref.parent / "test035b" / "test035b_gamma_tle" + if not dir_gamma.exists(): + dir_gamma = paths.output.parent / "test035b_gamma_tle" + + if dir_e.exists() and dir_gamma.exists(): + with tempfile.TemporaryDirectory() as tmpdir: + out_merge = pathlib.Path(tmpdir) / "merged" + res_merge = runner.invoke( + go, + [ + "--merge-only", + str(dir_e), + str(dir_gamma), + "-o", + str(out_merge), + "--e-factor", + "1.0", + ], + ) + print(f" --merge-only exit code: {res_merge.exit_code}") + merged_dose = out_merge / "output_dose.mhd" + merged_unc = out_merge / "output_dose_uncertainty.mhd" + if ( + res_merge.exit_code != 0 + or not merged_dose.exists() + or not merged_unc.exists() + ): + print(f" Merge failed or output missing: {res_merge.output}") + is_ok = False + else: + print(" Merge completed successfully, output images created.") + else: + print(f" Skipping --merge-only (test035b dirs not found: {dir_e})") + + print("\n=== Test 4: Short end-to-end analog CLI run with generic ion ===") + with tempfile.TemporaryDirectory() as tmpdir: + tmp = pathlib.Path(tmpdir) + cfg = { + "ct_image": str(dr_data / "29_CT_5mm_crop.mhd"), + "table_mat": str(dr_data / "Schneider2000MaterialsTable.txt"), + "table_density": str(dr_data / "Schneider2000DensitiesTable.txt"), + "activity_image": str(dr_data / "activity_test_crop_4mm.mhd"), + "density_tolerance_gcm3": 0.2, + "verbose": False, + } + cfg_path = tmp / "param.json" + cfg_path.write_text(json.dumps(cfg)) + out_sim = tmp / "output_cli" + + # Run with '89 225' (Ac225), low activity (10 Bq), 1 thread + res_sim = runner.invoke( + go, + [ + str(cfg_path), + "-r", + "89 225", + "-a", + "10", + "-t", + "1", + "-o", + str(out_sim), + ], + ) + print(f" CLI simulation exit code: {res_sim.exit_code}") + out_dose = out_sim / "output_dose.mhd" + out_stats = out_sim / "stats.txt" + if res_sim.exit_code != 0 or not out_dose.exists() or not out_stats.exists(): + print(f" Simulation failed or files missing! Output:\n{res_sim.output}") + is_ok = False + else: + print(f" Simulation produced: {[f.name for f in out_sim.glob('*.mhd')]}") + + utility.test_ok(is_ok) From 2eef1453b54e0a1fa876e472d2d5298c7cfed9e5 Mon Sep 17 00:00:00 2001 From: David Date: Thu, 10 Sep 2026 08:38:30 +0200 Subject: [PATCH 5/6] fix typo --- opengate/tests/src/source/test035a_dose_rate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/tests/src/source/test035a_dose_rate.py b/opengate/tests/src/source/test035a_dose_rate.py index b7eac05850..89c03c1185 100755 --- a/opengate/tests/src/source/test035a_dose_rate.py +++ b/opengate/tests/src/source/test035a_dose_rate.py @@ -60,7 +60,7 @@ print(h) is_ok = ( utility.assert_images( - paths.output_ref / " .mhd", + paths.output_ref / "edep.mhd", h.edep.get_output_path(), stats, tolerance=15, From a92a62bf132355ed95403e254107ee59d8511acf Mon Sep 17 00:00:00 2001 From: David Date: Thu, 10 Sep 2026 11:20:31 +0200 Subject: [PATCH 6/6] fix thread race issue --- .../opengate_lib/GateMaterialMuHandler.cpp | 60 +++++++------ .../opengate_lib/GateMaterialMuHandler.h | 2 - .../opengate_lib/GateMuTables.cpp | 84 ++++++++----------- .../opengate_core/opengate_lib/GateMuTables.h | 12 +-- .../opengate_lib/GateTLEDoseActor.cpp | 1 + .../test009_voxels_dynamic_jobs_merge.py | 0 6 files changed, 75 insertions(+), 84 deletions(-) mode change 100644 => 100755 opengate/tests/src/geometry/test009_voxels_dynamic_jobs_merge.py diff --git a/core/opengate_core/opengate_lib/GateMaterialMuHandler.cpp b/core/opengate_core/opengate_lib/GateMaterialMuHandler.cpp index 01fd96dcd7..8ac2588d23 100644 --- a/core/opengate_core/opengate_lib/GateMaterialMuHandler.cpp +++ b/core/opengate_core/opengate_lib/GateMaterialMuHandler.cpp @@ -9,6 +9,7 @@ #include "GateHelpers.h" #include "GateMuDatabase.h" #include "GateMuTables.h" +#include #include #include #include @@ -16,14 +17,15 @@ #include #include -// GateMaterialMuHandler *GateMaterialMuHandler::fSingletonMaterialMuHandler = -// nullptr; +G4Mutex sMaterialMuHandlerMutex = G4MUTEX_INITIALIZER; + std::map, std::shared_ptr> GateMaterialMuHandler::fInstances; std::shared_ptr GateMaterialMuHandler::GetInstance(std::string database, double energy_max) { + G4AutoLock mutex(&sMaterialMuHandlerMutex); // Create a key based on database and energy_max auto key = std::make_tuple(database, energy_max); @@ -49,55 +51,64 @@ GateMaterialMuHandler::GateMaterialMuHandler() { fEnergyNumber = 40; fAtomicShellEnergyMin = 1. * CLHEP::keV; fPrecision = 0.01; - fLastCouple = nullptr; - fLastMuTable = nullptr; } GateMaterialMuHandler::~GateMaterialMuHandler() { delete[] fElementsTable; } -void GateMaterialMuHandler::CheckLastCall(const G4MaterialCutsCouple *couple) { +void GateMaterialMuHandler::CheckLastCall(const G4MaterialCutsCouple *) { if (!fIsInitialized) { Initialize(); } - if (couple != fLastCouple) { - fLastCouple = couple; - fLastMuTable = fCoupleTable[fLastCouple]; +} + +GateMuTable * +GateMaterialMuHandler::GetMuTable(const G4MaterialCutsCouple *couple) { + if (!fIsInitialized) { + Initialize(); } + thread_local const GateMaterialMuHandler *lastHandler = nullptr; + thread_local const G4MaterialCutsCouple *lastCouple = nullptr; + thread_local GateMuTable *lastTable = nullptr; + if (this == lastHandler && couple == lastCouple && lastTable != nullptr) { + return lastTable; + } + auto it = fCoupleTable.find(couple); + if (it != fCoupleTable.end()) { + lastHandler = this; + lastCouple = couple; + lastTable = it->second; + return lastTable; + } + return nullptr; } double GateMaterialMuHandler::GetDensity(const G4MaterialCutsCouple *couple) { - CheckLastCall(couple); - return fLastMuTable->GetDensity(); + auto *muTable = GetMuTable(couple); + return muTable ? muTable->GetDensity() : 0.0; } double GateMaterialMuHandler::GetMuEnOverRho(const G4MaterialCutsCouple *couple, double energy) { - CheckLastCall(couple); - return fLastMuTable->GetMuEnOverRho(energy); + auto *muTable = GetMuTable(couple); + return muTable ? muTable->GetMuEnOverRho(energy) : 0.0; } double GateMaterialMuHandler::GetMuEn(const G4MaterialCutsCouple *couple, double energy) { - CheckLastCall(couple); - return fLastMuTable->GetMuEn(energy); + auto *muTable = GetMuTable(couple); + return muTable ? muTable->GetMuEn(energy) : 0.0; } double GateMaterialMuHandler::GetMuOverRho(const G4MaterialCutsCouple *couple, double energy) { - CheckLastCall(couple); - return fLastMuTable->GetMuOverRho(energy); + auto *muTable = GetMuTable(couple); + return muTable ? muTable->GetMuOverRho(energy) : 0.0; } double GateMaterialMuHandler::GetMu(const G4MaterialCutsCouple *couple, double energy) { - CheckLastCall(couple); - return fLastMuTable->GetMu(energy); -} - -GateMuTable * -GateMaterialMuHandler::GetMuTable(const G4MaterialCutsCouple *couple) { - CheckLastCall(couple); - return fLastMuTable; + auto *muTable = GetMuTable(couple); + return muTable ? muTable->GetMu(energy) : 0.0; } inline double interpolation(double Xa, double Xb, double Ya, double Yb, @@ -106,6 +117,7 @@ inline double interpolation(double Xa, double Xb, double Ya, double Yb, } void GateMaterialMuHandler::Initialize() { + G4AutoLock mutex(&sMaterialMuHandlerMutex); if (fIsInitialized) return; diff --git a/core/opengate_core/opengate_lib/GateMaterialMuHandler.h b/core/opengate_core/opengate_lib/GateMaterialMuHandler.h index f939584182..a8afad7255 100644 --- a/core/opengate_core/opengate_lib/GateMaterialMuHandler.h +++ b/core/opengate_core/opengate_lib/GateMaterialMuHandler.h @@ -109,8 +109,6 @@ class GateMaterialMuHandler { int fEnergyNumber; double fAtomicShellEnergyMin; double fPrecision; - const G4MaterialCutsCouple *fLastCouple; - GateMuTable *fLastMuTable; }; #endif diff --git a/core/opengate_core/opengate_lib/GateMuTables.cpp b/core/opengate_core/opengate_lib/GateMuTables.cpp index fd156dcc2e..ec245d0098 100644 --- a/core/opengate_core/opengate_lib/GateMuTables.cpp +++ b/core/opengate_core/opengate_lib/GateMuTables.cpp @@ -12,10 +12,6 @@ GateMuTable::GateMuTable(const G4MaterialCutsCouple *couple, G4int size) { fMu = new double[size]; fMuEn = new double[size]; fSize = size; - fLastMu = -1.0; - fLastMuEn = -1.0; - fLastEnergyMu = -1.0; - fLastEnergyMuEn = -1.0; mCouple = couple; mDensity = -1; @@ -44,69 +40,57 @@ inline double interpol(double x1, double x, double x2, double y1, double y2) { // storage } -double GateMuTable::GetMuEnOverRho(double energy) { - if (energy != fLastEnergyMuEn) { - fLastEnergyMuEn = energy; - - energy = log(energy); - - int inf = 0; - int sup = fSize - 1; - while (sup - inf > 1) { - int tmp_bound = (inf + sup) / 2; - if (fEnergy[tmp_bound] > energy) { - sup = tmp_bound; - } else { - inf = tmp_bound; - } - } - double e_inf = fEnergy[inf]; - double e_sup = fEnergy[sup]; +double GateMuTable::GetMuEnOverRho(double energy) const { + const double log_energy = log(energy); - if (energy > e_inf && energy < e_sup) { - fLastMuEn = exp(interpol(e_inf, energy, e_sup, fMuEn[inf], fMuEn[sup])); + int inf = 0; + int sup = fSize - 1; + while (sup - inf > 1) { + int tmp_bound = (inf + sup) / 2; + if (fEnergy[tmp_bound] > log_energy) { + sup = tmp_bound; } else { - fLastMuEn = exp(fMuEn[inf]); + inf = tmp_bound; } } + double e_inf = fEnergy[inf]; + double e_sup = fEnergy[sup]; - return fLastMuEn; + if (log_energy > e_inf && log_energy < e_sup) { + return exp(interpol(e_inf, log_energy, e_sup, fMuEn[inf], fMuEn[sup])); + } else { + return exp(fMuEn[inf]); + } } -double GateMuTable::GetMuEn(double energy) { +double GateMuTable::GetMuEn(double energy) const { return (GetMuEnOverRho(energy) * mDensity); } -double GateMuTable::GetMuOverRho(double energy) { - if (energy != fLastEnergyMu) { - fLastEnergyMu = energy; - - energy = log(energy); - - int inf = 0; - int sup = fSize - 1; - while (sup - inf > 1) { - int tmp_bound = (inf + sup) / 2; - if (fEnergy[tmp_bound] > energy) { - sup = tmp_bound; - } else { - inf = tmp_bound; - } - } - double e_inf = fEnergy[inf]; - double e_sup = fEnergy[sup]; +double GateMuTable::GetMuOverRho(double energy) const { + const double log_energy = log(energy); - if (energy > e_inf && energy < e_sup) { - fLastMu = exp(interpol(e_inf, energy, e_sup, fMu[inf], fMu[sup])); + int inf = 0; + int sup = fSize - 1; + while (sup - inf > 1) { + int tmp_bound = (inf + sup) / 2; + if (fEnergy[tmp_bound] > log_energy) { + sup = tmp_bound; } else { - fLastMu = exp(fMu[inf]); + inf = tmp_bound; } } + double e_inf = fEnergy[inf]; + double e_sup = fEnergy[sup]; - return fLastMu; + if (log_energy > e_inf && log_energy < e_sup) { + return exp(interpol(e_inf, log_energy, e_sup, fMu[inf], fMu[sup])); + } else { + return exp(fMu[inf]); + } } -double GateMuTable::GetMu(double energy) { +double GateMuTable::GetMu(double energy) const { return (GetMuOverRho(energy) * mDensity); } diff --git a/core/opengate_core/opengate_lib/GateMuTables.h b/core/opengate_core/opengate_lib/GateMuTables.h index 176de85cc3..90f1492ac5 100644 --- a/core/opengate_core/opengate_lib/GateMuTables.h +++ b/core/opengate_core/opengate_lib/GateMuTables.h @@ -20,13 +20,13 @@ class GateMuTable { void PutValue(int index, double energy, double mu, double mu_en) const; - double GetMuEn(double energy); + double GetMuEn(double energy) const; - double GetMuEnOverRho(double energy); + double GetMuEnOverRho(double energy) const; - double GetMu(double energy); + double GetMu(double energy) const; - double GetMuOverRho(double energy); + double GetMuOverRho(double energy) const; const G4MaterialCutsCouple *GetMaterialCutsCouple() const; @@ -48,10 +48,6 @@ class GateMuTable { double *fEnergy; double *fMu; double *fMuEn; - double fLastEnergyMu; - double fLastEnergyMuEn; - double fLastMu; - double fLastMuEn; G4int fSize; }; diff --git a/core/opengate_core/opengate_lib/GateTLEDoseActor.cpp b/core/opengate_core/opengate_lib/GateTLEDoseActor.cpp index f6566dd795..29f0fd14b8 100644 --- a/core/opengate_core/opengate_lib/GateTLEDoseActor.cpp +++ b/core/opengate_core/opengate_lib/GateTLEDoseActor.cpp @@ -129,6 +129,7 @@ void GateTLEDoseActor::BeginOfEventAction(const G4Event *event) { fMaterialMuHandler = GateMaterialMuHandler::GetInstance(fDatabase, 5 * CLHEP::MeV); } + fMaterialMuHandler->Initialize(); } } diff --git a/opengate/tests/src/geometry/test009_voxels_dynamic_jobs_merge.py b/opengate/tests/src/geometry/test009_voxels_dynamic_jobs_merge.py old mode 100644 new mode 100755