From b759a7f7de877aa3374e9a7e4f0024734396238b Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 7 Nov 2023 15:49:14 +0100 Subject: [PATCH 01/37] correct arf output size (was wrong) --- opengate/actors/arfactors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/actors/arfactors.py b/opengate/actors/arfactors.py index 8dd9d9613..b1f62b03f 100644 --- a/opengate/actors/arfactors.py +++ b/opengate/actors/arfactors.py @@ -268,7 +268,7 @@ def EndSimulationAction(self): # Should we keep the first slice (with all hits) ? if not self.user_info.enable_hit_slice: self.output_image = self.output_image[1:, :, :] - self.param.image_size[1] = self.param.image_size[1] - 1 + self.param.image_size[0] = self.param.image_size[0] - 1 # convert to itk image self.output_image = itk.image_from_array(self.output_image) From d9378747b72cb46d191ff75bcf8b4ab33fc36305 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 7 Nov 2023 15:50:44 +0100 Subject: [PATCH 02/37] adapt to new API --- opengate/contrib/spect/genm670.py | 5 ++--- opengate/tests/src/test028_ge_nm670_spect_1_colli_visu.py | 3 +-- opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py | 7 ++++--- opengate/tests/src/test045_speedup.py | 2 +- 4 files changed, 8 insertions(+), 9 deletions(-) diff --git a/opengate/contrib/spect/genm670.py b/opengate/contrib/spect/genm670.py index 90cfcce60..df7f2e902 100644 --- a/opengate/contrib/spect/genm670.py +++ b/opengate/contrib/spect/genm670.py @@ -266,8 +266,6 @@ def hegp_collimator_repeater(sim, name, core, debug): holep = RepeatParametrisedVolume(repeated_volume=hole) holep.linear_repeat = [54, 70, 1] if debug: - holep.linear_repeat = [54, 70, 1] - else: holep.linear_repeat = [10, 10, 1] holep.translation = [10.0459 * mm, 5.8 * mm, 0] # dot it twice, with the following offset @@ -307,12 +305,13 @@ def lehr_collimator_repeater(sim, name, core, debug): hole.radius = 0.075 * cm hole.material = "G4_AIR" hole.mother = core.name + hole.build_physical_volume = False sim.volume_manager.add_volume(hole) # parameterised holes holep = RepeatParametrisedVolume(repeated_volume=hole) if debug: - holep.linear_repeat = [10, 10, 1] + holep.linear_repeat = [30, 30, 1] else: holep.linear_repeat = [183, 235, 1] diff --git a/opengate/tests/src/test028_ge_nm670_spect_1_colli_visu.py b/opengate/tests/src/test028_ge_nm670_spect_1_colli_visu.py index 56b13c195..4e2141140 100755 --- a/opengate/tests/src/test028_ge_nm670_spect_1_colli_visu.py +++ b/opengate/tests/src/test028_ge_nm670_spect_1_colli_visu.py @@ -15,8 +15,7 @@ ui = sim.user_info ui.g4_verbose = False ui.visu = True - ui.visu_type = "vrml" - ui.visu_filename = "geant4VisuFile.wrl" + ui.visu_type = "qt" ui.number_of_threads = 1 ui.check_volumes_overlap = False diff --git a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py index e9858a5ff..a7e8609d8 100644 --- a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py +++ b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py @@ -67,7 +67,7 @@ def create_pet_simulation(sim, param): ) # physic list - sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4" + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" sim.physics_manager.enable_decay = False # p.apply_cuts = True @@ -93,7 +93,7 @@ def add_analytical_phantom(sim, param): def add_voxelized_phantom(sim, param): print("Phantom: IEC voxelized: ", param.iec_vox_mhd) iec = sim.add_volume("Image", "iec") - gate_iec.create_material() + gate_iec.create_material(sim) iec.image = param.iec_vox_mhd iec.material = "G4_AIR" labels = json.loads(open(param.iec_vox_json).read()) @@ -131,6 +131,7 @@ def add_pet(sim, param): # hits collection hc = sim.add_actor("DigitizerHitsCollectionActor", "Hits") + crystal = None # get crystal volume by looking for the word crystal in the name for k, v in sim.volume_manager.volumes.items(): if "crystal" in k: @@ -331,7 +332,7 @@ def add_voxelized_source(sim, p): # compute volume to convert Bqml in Bq img = itk.imread(p.source_vox_mhd) stats = Box(gt.imageStatistics(img, None, False, 10)) - info = gate.get_info_from_image(img) + info = gate.image.get_info_from_image(img) vol = stats.sum * info.spacing[0] * info.spacing[1] * info.spacing[2] * mm3 print(f"Volume source {vol} mm3") ac = p.activity_Bqml * vol / cm3 diff --git a/opengate/tests/src/test045_speedup.py b/opengate/tests/src/test045_speedup.py index 898a5f48f..3778eeb10 100755 --- a/opengate/tests/src/test045_speedup.py +++ b/opengate/tests/src/test045_speedup.py @@ -63,7 +63,7 @@ def run_test_045_speedrun( p.source_vox_mhd = str(paths.data / "iec_source_4mm.mhd") p.gaga_pth = paths.data / "pth120_test9221_GP_0GP_10.0_100000.pth" - gate.helpers.print_dic(p) + gate.utility.print_dic(p) # output if output_folder == "AUTO": From 0f2c2e17fc64df817f02501eb2636a79df3cf08b Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 7 Nov 2023 15:51:59 +0100 Subject: [PATCH 03/37] option to use image origin --- opengate/sources/gansources.py | 35 +++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/opengate/sources/gansources.py b/opengate/sources/gansources.py index 2e59ead35..3e85f9bda 100644 --- a/opengate/sources/gansources.py +++ b/opengate/sources/gansources.py @@ -188,23 +188,41 @@ def sample_indices_phys(self, n, rs=np.random): class VoxelizedSourceConditionGenerator: - def __init__(self, activity_source_filename, rs=np.random): + def __init__( + self, activity_source_filename, rs=np.random, use_activity_origin=False + ): self.activity_source_filename = str(activity_source_filename) + # options + self.compute_directions = False + self.use_activity_origin = use_activity_origin + # variables self.image = None self.cdf_x = self.cdf_y = self.cdf_z = None self.rs = rs self.img_info = None self.sampler = None - self.initialize_source() - self.compute_directions = False + self.points_offset = None + # init + self.is_initialized = False def initialize_source(self): self.image = itk.imread(self.activity_source_filename) self.img_info = get_info_from_image(self.image) self.sampler = VoxelizedSourcePDFSampler(self.image) self.rs = np.random + # we set the points in the g4 coord system (according to the center of the image) + # or according to the activity source image origin + if self.use_activity_origin is True: + self.points_offset = self.img_info.origin + else: + hs = self.img_info.spacing / 2.0 + self.points_offset = -hs * self.img_info.size + hs + self.is_initialized = True def generate_condition(self, n): + if self.is_initialized is False: + self.initialize_source() + # i j k is in np array order = z y x # but img_info is in the order x y z i, j, k = self.sampler.sample_indices(n, self.rs) @@ -223,15 +241,14 @@ def generate_condition(self, n): z = self.img_info.spacing[0] * i + rx # x,y,z are in the image coord system - # we set in the g4 coord system: according to the center of the image - p = np.column_stack((x, y, z)) - hs * self.img_info.size + hs + # tey are offset according to the coord system (image center or image offset) + p = np.column_stack((x, y, z)) + self.points_offset # need direction ? - if self.compute_directions: - v = generate_isotropic_directions(n, rs=self.rs) - return np.column_stack((p, v)) - else: + if self.compute_directions is False: return p + v = generate_isotropic_directions(n, rs=self.rs) + return np.column_stack((p, v)) class GANSource(GenericSource): From 23cc69f143a78a605474589de7eac6f68860f252 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 13 Nov 2023 09:48:10 +0100 Subject: [PATCH 04/37] Add test 066 --- .../test066_spect_gaga_garf_1_reference.py | 56 +++++ .../tests/src/test066_spect_gaga_garf_2.py | 73 ++++++ .../test066_spect_gaga_garf_3_standalone.py | 79 ++++++ .../src/test066_spect_gaga_garf_4_analyse1.py | 60 +++++ .../src/test066_spect_gaga_garf_5_analyse2.py | 60 +++++ .../src/test066_spect_gaga_garf_helpers.py | 238 ++++++++++++++++++ 6 files changed, 566 insertions(+) create mode 100755 opengate/tests/src/test066_spect_gaga_garf_1_reference.py create mode 100755 opengate/tests/src/test066_spect_gaga_garf_2.py create mode 100755 opengate/tests/src/test066_spect_gaga_garf_3_standalone.py create mode 100755 opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py create mode 100755 opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py create mode 100644 opengate/tests/src/test066_spect_gaga_garf_helpers.py diff --git a/opengate/tests/src/test066_spect_gaga_garf_1_reference.py b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py new file mode 100755 index 000000000..7f2452169 --- /dev/null +++ b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from test066_spect_gaga_garf_helpers import * +from opengate.tests import utility + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, None, "test066") + output_path = paths.output + + # create the simulation + sim = gate.Simulation() + simu_name = "test066_1_reference" + + # options + ui = sim.user_info + ui.number_of_threads = 10 + # ui.visu = True + ui.visu_type = "vrml" + ui.random_seed = "auto" + ui.check_volumes_overlap = False + + # units + mm = gate.g4_units.mm + sec = gate.g4_units.second + Bq = gate.g4_units.Bq + cm3 = gate.g4_units.cm3 + BqmL = Bq / cm3 + + # main elements : spect + phantom + proj1, proj2 = create_simu_with_genm670(sim, debug=ui.visu) + proj1.output = f"{output_path}/{simu_name}_0.mhd" + proj2.output = f"{output_path}/{simu_name}_1.mhd" + + # add IEC phantom + gate_iec.add_iec_phantom(sim, name="iec") + sim.set_production_cut("iec", "all", 1 * mm) + + # sources IEC + ac = 1e5 * BqmL # 1e5 = about 10 min with 10 threads linux + if ui.visu: + ac = 0.01 * BqmL + total_activity, w, e = add_iec_Tc99m_source(sim, ac) + + # duration + set_duration(sim, total_activity, w, 30 * sec) + + # run + sim.run() + + # print results at the end + stats = sim.output.get_actor("stats") + stats.write(f"{output_path}/{simu_name}_stats.txt") + print(stats) + print(output_path) diff --git a/opengate/tests/src/test066_spect_gaga_garf_2.py b/opengate/tests/src/test066_spect_gaga_garf_2.py new file mode 100755 index 000000000..abf22d43e --- /dev/null +++ b/opengate/tests/src/test066_spect_gaga_garf_2.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from test066_spect_gaga_garf_helpers import * +from opengate.tests import utility + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, None, "test066") + output_path = paths.output + + # create the simulation + sim = gate.Simulation() + simu_name = "test066_2" + + # options + ui = sim.user_info + ui.number_of_threads = 5 + # ui.visu = True + ui.visu_type = "vrml" + ui.random_seed = "auto" + # ui.check_volumes_overlap = False + + # units + mm = gate.g4_units.mm + sec = gate.g4_units.second + Bq = gate.g4_units.Bq + cm3 = gate.g4_units.cm3 + BqmL = Bq / cm3 + + # activity + # FIXME to compute exact number + w, e = get_rad_gamma_energy_spectrum("Tc99m") + # Estimated gammas 127,008,708 gammas (weights = 0.8850) + # so, because we use ARF, about 1/2 particles needed + total_activity = 127008708 / 30 * Bq / ui.number_of_threads / 2 + print(f"Total activity: {total_activity/Bq:.0f} Bq") + if ui.visu: + total_activity = 10 * Bq + + # source + # voxelize_iec_phantom -o ../data/iec_2mm.mhd -s 2 --output_source data/iec_source_same_spheres_2mm.mhd + # -a 1 1 1 1 1 1 --no_shell + + # main elements : spect + phantom + activity_source = paths.data / "iec_source_same_spheres_2mm.mhd" + gaga_pth_filename = ( + paths.data + / "gate" + / "gate_test038_gan_phsp_spect" + / "pth2" + / "test001_GP_0GP_10_50000.pth" + ) + garf_pth_filename = ( + paths.data / "gate" / "gate_test043_garf" / "data" / "pth" / "arf_Tc99m_v3.pth" + ) + arf1, arf2 = create_simu_with_gaga( + sim, total_activity, activity_source, gaga_pth_filename, garf_pth_filename + ) + arf1.output = f"{output_path}/{simu_name}_0.mhd" + arf2.output = f"{output_path}/{simu_name}_1.mhd" + + # duration + set_duration(sim, total_activity, w, 30 * sec) + + # run + sim.run() + + # print results at the end + stats = sim.output.get_actor("stats") + stats.write(f"{output_path}/{simu_name}_stats.txt") + print(stats) + print(output_path) diff --git a/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py b/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py new file mode 100755 index 000000000..44e540f92 --- /dev/null +++ b/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from opengate.tests import utility +import gaga_phsp as gaga +import itk +from box import Box +from scipy.spatial.transform import Rotation +import opengate as gate +from opengate.contrib.spect import genm670 + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, None, "test066") + output_path = paths.output + + # create the simulation + simu_name = "test066_3_standalone" + g = paths.data / "gate" + + # info about ge nm670 + mm = gate.g4_units.mm + head_radius = 280 * mm + pos, crystal_dist, psd = genm670.compute_plane_position_and_distance_to_crystal( + "lehr" + ) + + # parameters + activity_source = paths.data / "iec_source_same_spheres_2mm.mhd" + gaga_pth_filename = ( + g / "gate_test038_gan_phsp_spect" / "pth2" / "test001_GP_0GP_10_50000.pth" + ) + garf_pth_filename = g / "gate_test043_garf" / "data" / "pth" / "arf_Tc99m_v3.pth" + gaga_user_info = Box( + { + "pth_filename": gaga_pth_filename, + "activity_source": activity_source, + "batch_size": 1e5, + "gpu_mode": "auto", + "backward_distance": 50 * mm, + } + ) + garf_user_info = Box( + { + "pth_filename": garf_pth_filename, + "image_size": [128, 128], + "image_spacing": [3 * mm, 3 * mm], + "plane_distance": -head_radius - pos, + "distance_to_crystal": crystal_dist, + "batch_size": 1e5, + "gpu_mode": "auto", + } + ) + + print(f"{garf_user_info.plane_distance=}") + + # initialize gaga and garf (read the NN) + gaga.gaga_garf_generate_spect_initialize(gaga_user_info, garf_user_info) + + # Initial rotation angle of the head (identity here) + r = Rotation.from_euler("y", 0, degrees=True) + garf_user_info.plane_rotation = r + + # GO + n = 127008708 / 20 + angle_rotations = [ + Rotation.from_euler("y", 0, degrees=True), + Rotation.from_euler("y", 180, degrees=True), + ] + images = gaga.gaga_garf_generate_spect( + gaga_user_info, garf_user_info, n, angle_rotations + ) + + # save image + i = 0 + for image in images: + output = output_path / f"{simu_name}_{i}.mhd" + print(f"Done, saving output in {output}") + itk.imwrite(image, output) + i += 1 diff --git a/opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py b/opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py new file mode 100755 index 000000000..813c0a2f2 --- /dev/null +++ b/opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from opengate.tests import utility +from box import Box + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, None, "test066") + output_path = paths.output + + # folders and names + ref_simu_name = "test066_1_reference" + ref_output_folder = output_path + test_simu_name = "test066_2" + test_output_folder = output_path + + # image filenames + ref_names = [ + ref_output_folder / f"{ref_simu_name}_0.mhd", + ref_output_folder / f"{ref_simu_name}_1.mhd", + ] + test_names = [ + test_output_folder / f"{test_simu_name}_0.mhd", + test_output_folder / f"{test_simu_name}_1.mhd", + ] + + options = Box() + options.scaling = 2 + options.n_slice = 1 + options.window_width = 33 + options.window_level = 17 + options.crop_center = (63, 73) + options.crop_width = (100, 100) + options.hline = 58 + options.vline = 44 + options.width = 10 + options.lab_ref = "GAGA-GARF in Gate" + options.lab_test = "GAGA-GARF standalone" + options.title = "Ref vs GAGA-REF in Gate" + + plt = utility.plot_compare_profile(ref_names, test_names, options) + + f = test_output_folder / f"{test_simu_name}.pdf" + print("Save figure in ", f) + plt.savefig(f, bbox_inches="tight", format="pdf") + + # test + print() + is_ok = True + for ref_name, test_name in zip(ref_names, test_names): + is_ok = utility.assert_images( + ref_name, + test_name, + axis="x", + scaleImageValuesFactor=options.scaling, + sum_tolerance=8, + tolerance=110, + ) + + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py b/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py new file mode 100755 index 000000000..0c77182ae --- /dev/null +++ b/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from opengate.tests import utility +from box import Box + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, None, "test066") + output_path = paths.output + + # folders and names + ref_simu_name = "test066_2" + ref_output_folder = output_path + test_simu_name = "test066_3_standalone" + test_output_folder = output_path + + # image filenames + ref_names = [ + ref_output_folder / f"{ref_simu_name}_0.mhd", + ref_output_folder / f"{ref_simu_name}_1.mhd", + ] + test_names = [ + test_output_folder / f"{test_simu_name}_0.mhd", + test_output_folder / f"{test_simu_name}_1.mhd", + ] + + options = Box() + options.scaling = 10 + options.n_slice = 1 + options.window_width = 33 + options.window_level = 17 + options.crop_center = (63, 73) + options.crop_width = (100, 100) + options.hline = 58 + options.vline = 44 + options.width = 10 + options.lab_ref = "GAGA-GARF in Gate" + options.lab_test = "GAGA-GARF standalone" + options.title = "GAGA-REF in Gate vs standalone" + + plt = utility.plot_compare_profile(ref_names, test_names, options) + + f = test_output_folder / f"{test_simu_name}.pdf" + print("Save figure in ", f) + plt.savefig(f, bbox_inches="tight", format="pdf") + + # test + print() + is_ok = True + for ref_name, test_name in zip(ref_names, test_names): + is_ok = utility.assert_images( + ref_name, + test_name, + axis="x", + scaleImageValuesFactor=options.scaling, + sum_tolerance=15, + tolerance=105, + ) + + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test066_spect_gaga_garf_helpers.py b/opengate/tests/src/test066_spect_gaga_garf_helpers.py new file mode 100644 index 000000000..ffe83a3b3 --- /dev/null +++ b/opengate/tests/src/test066_spect_gaga_garf_helpers.py @@ -0,0 +1,238 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import opengate.contrib.phantoms.nemaiec as gate_iec +import opengate.sources.gansources as gansources +from opengate.sources.generic import get_rad_gamma_energy_spectrum +from opengate.contrib.spect import genm670 +from scipy.spatial.transform import Rotation + + +def create_world(sim): + # world size + m = gate.g4_units.m + sim.world.size = [1 * m, 1 * m, 1 * m] + sim.world.material = "G4_AIR" + + +def set_phys(sim): + m = gate.g4_units.m + p = sim.get_physics_user_info() + p.physics_list_name = "G4EmStandardPhysics_option3" + sim.set_production_cut("world", "all", 1e3 * m) + + +def create_simu_with_genm670(sim, collimator_type="lehr", debug=False): + # units + mm = gate.g4_units.mm + + # world size + create_world(sim) + + # spect system + head1, crystal1 = genm670.add_ge_nm67_spect_head( + sim, "spect1", collimator_type=collimator_type, debug=debug + ) + head2, crystal2 = genm670.add_ge_nm67_spect_head( + sim, "spect2", collimator_type=collimator_type, debug=debug + ) + + # default rotation to be in front of the IEC (not a realistic orientation) + head1.translation = [0, 0, -280 * mm] + head2.translation = [0, 0, 280 * mm] + head1.rotation = Rotation.from_euler("y", 0, degrees=True).as_matrix() + head2.rotation = Rotation.from_euler("y", 180, degrees=True).as_matrix() + + # digitizer + keV = gate.g4_units.keV + channels = [ + {"name": f"scatter", "min": 114 * keV, "max": 126 * keV}, + {"name": f"peak140", "min": 126 * keV, "max": 154 * keV}, + ] + proj1 = genm670.add_digitizer(sim, crystal1.name, channels) + proj2 = genm670.add_digitizer(sim, crystal2.name, channels) + proj1.spacing = [3 * mm, 3 * mm] + proj2.spacing = [3 * mm, 3 * mm] + + # stats + sim.add_actor("SimulationStatisticsActor", "stats") + + # phys + set_phys(sim) + sim.set_production_cut("spect1", "all", 1 * mm) + sim.set_production_cut("spect2", "all", 1 * mm) + + return proj1, proj2 + + +def create_simu_with_gaga( + sim, total_activity, activity_source, gaga_pth_filename, garf_pth_filename +): + # world size + create_world(sim) + + # add gaga source + add_gaga_source(sim, total_activity, activity_source, gaga_pth_filename) + + # add arf plane + mm = gate.g4_units.mm + det1, det2 = add_garf_detector_planes(sim, head_radius=280 * mm) + arf1, arf2 = add_garf_detectors(sim, garf_pth_filename, det1, det2) + + # stats + sim.add_actor("SimulationStatisticsActor", "stats") + + # phys + set_phys(sim) + + return arf1, arf2 + + +def add_gaga_source(sim, total_activity, activity_source, gaga_pth_filename): + mm = gate.g4_units.mm + keV = gate.g4_units.keV + gsource = sim.add_source("GANSource", "gaga") + gsource.particle = "gamma" + gsource.activity = total_activity + gsource.pth_filename = gaga_pth_filename + gsource.position_keys = ["PrePosition_X", "PrePosition_Y", "PrePosition_Z"] + gsource.backward_distance = 50 * mm + gsource.backward_force = True + gsource.direction_keys = ["PreDirection_X", "PreDirection_Y", "PreDirection_Z"] + gsource.energy_key = "KineticEnergy" + gsource.energy_min_threshold = 10 * keV + gsource.skip_policy = "ZeroEnergy" + gsource.weight_key = None + gsource.time_key = None + gsource.batch_size = 5e4 # OSX + # Linux 5e4 with 1 thread, 1e4 with 8 threads + gsource.batch_size = 2e5 / sim.user_info.number_of_threads + gsource.verbose_generator = False # True + gsource.gpu_mode = "auto" + + # conditional generator + cond_generator = gansources.VoxelizedSourceConditionGenerator( + activity_source, use_activity_origin=True + ) + cond_generator.compute_directions = True + gene = gansources.GANSourceConditionalGenerator( + gsource, cond_generator.generate_condition + ) + gsource.generator = gene + + +def add_garf_detector_planes(sim, head_radius): + cm = gate.g4_units.cm + nm = gate.g4_units.nm + # GE colli size + colli_size = [54.6 * cm, 40.6 * cm] # FIXME function + # colli_size = [100 * cm, 100 * cm] # FIXME function + colli_size = [128 * 3, 128 * 3] # FIXME function + pos, crystal_dist, psd = genm670.compute_plane_position_and_distance_to_crystal( + "lehr" + ) + print(pos, crystal_dist, psd) + print(f"The detector plane position = {pos} mm") + print(f"The detector head radius = {head_radius} mm") + + detector_plane1 = sim.add_volume("Box", "det_plane1") + detector_plane1.material = "G4_Galactic" + detector_plane1.color = [1, 0, 0, 1] + detector_plane1.size = [colli_size[0], colli_size[1], 1 * nm] + r = Rotation.from_euler("x", 0, degrees=True) + detector_plane1.rotation = r.as_matrix() + detector_plane1.translation = [0, 0, -head_radius - pos] + print(f"{detector_plane1.translation=}") + + detector_plane2 = sim.add_volume("Box", "det_plane2") + detector_plane2.material = "G4_Galactic" + detector_plane2.color = [1, 0, 0, 1] + detector_plane2.size = [colli_size[0], colli_size[1], 1 * nm] + r = Rotation.from_euler("y", 180, degrees=True) + detector_plane2.rotation = r.as_matrix() + detector_plane2.translation = [0, 0, -head_radius - pos] + detector_plane2.translation = r.apply(detector_plane2.translation) + print(f"{detector_plane2.translation=}") + # FIXME apply rotation to translation + + return detector_plane1, detector_plane2 + + +def add_garf_detectors(sim, pth_filename, det1, det2): + mm = gate.g4_units.mm + detectors = [det1, det2] + + pos, crystal_dist, psd = genm670.compute_plane_position_and_distance_to_crystal( + "lehr" + ) + print(f"{pos=}") + print(f"{crystal_dist=}") + print(f"{psd=}") + + arfs = [] + for det in detectors: + # arf actor + arf = sim.add_actor("ARFActor", f"arf_{det.name}") + arf.mother = det.name + arf.batch_size = 1e5 + arf.image_size = [128, 128] + arf.image_spacing = [3 * mm, 3 * mm] + arf.verbose_batch = True + arf.distance_to_crystal = crystal_dist + arf.gpu_mode = "auto" + arf.pth_filename = pth_filename + arf.enable_hit_slice = False + arfs.append(arf) + + return arfs[0], arfs[1] + + +def set_duration(sim, total_activity, w, duration): + Bq = gate.g4_units.Bq + sec = gate.g4_units.second + ui = sim.user_info + nb_decays = total_activity / Bq * duration / sec * ui.number_of_threads + weights = sum(w) + print(f"Estimated total decay {nb_decays:.0f} decays") + print( + f"Estimated gammas {nb_decays * weights:.0f} gammas (weights = {weights:.4f})" + ) + sim.run_timing_intervals = [[0, duration]] + + +def add_iec_Tc99m_source(sim, activity_concentration): + cm3 = gate.g4_units.cm3 + Bq = gate.g4_units.Bq + kBq = 1000 * Bq + MBq = 1000 * kBq + BqmL = Bq / cm3 + ui = sim.user_info + + # same concentration in all spheres + a = activity_concentration / ui.number_of_threads + all_activities = [a] * 6 + + # Tc99m + sources = gate_iec.add_spheres_sources(sim, "iec", "src", "all", all_activities) + total_activity = 0 + print(f"Activity concentration = {all_activities[5] / BqmL:.0f} BqmL") + w, e = get_rad_gamma_energy_spectrum("Tc99m") + for source in sources: + source.particle = "gamma" + source.energy.type = "spectrum_lines" + source.energy.spectrum_weight = w + source.energy.spectrum_energy = e + total_activity += source.activity + print( + f"Activity {source.name} = {source.activity / Bq:.2f} Bq {activity_concentration / BqmL:.0f} Bq/mL" + ) + + print(f"Total activity 1 thread = {total_activity / Bq:.2f} Bq") + print(f"Total activity 1 thread = {total_activity / kBq:.2f} kBq") + print(f"Total activity 1 thread = {total_activity / MBq:.2f} MBq") + print( + f"Total activity {ui.number_of_threads} threads = {total_activity / MBq * ui.number_of_threads:.2f} MBq" + ) + + return total_activity, w, e From 194fb5bb7712cca3fe70512cf96ba7b53bf22225 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 13 Nov 2023 09:49:26 +0100 Subject: [PATCH 05/37] add plot for test66 --- opengate/actors/arfactors.py | 26 +- opengate/bin/opengate_tests | 9 +- opengate/contrib/spect/genm670.py | 17 +- opengate/tests/data | 2 +- .../src/test045_gan_phsp_pet_gan_helpers.py | 4 +- opengate/tests/utility.py | 247 +++++++++++++++++- 6 files changed, 287 insertions(+), 18 deletions(-) diff --git a/opengate/actors/arfactors.py b/opengate/actors/arfactors.py index b1f62b03f..87d9391ec 100644 --- a/opengate/actors/arfactors.py +++ b/opengate/actors/arfactors.py @@ -115,6 +115,8 @@ def __init__(self, user_info): ActorBase.__init__(self, user_info) g4.GateARFActor.__init__(self, user_info.__dict__) # import module + self.debug_nb_hits_before = None + self.debug_nb_hits = 0 self.garf = import_garf() if self.garf is None: print("Cannot run GANSource") @@ -155,13 +157,17 @@ def initialize(self, volume_engine=None): self.ActorInitialize() self.SetARFFunction(self.apply) self.user_info.output_image = None + self.debug_nb_hits_before = 0 + self.debug_nb_hits = 0 # load the pth file self.nn, self.model = self.garf.load_nn(self.pth_filename, verbose=False) p = self.param p.batch_size = int(float(self.user_info.batch_size)) - # size and spacing (2D) + # size and spacing (2D) (force to float) + self.user_info.image_spacing[0] = float(self.user_info.image_spacing[0]) + self.user_info.image_spacing[1] = float(self.user_info.image_spacing[1]) p.image_size = self.user_info.image_size p.image_spacing = self.user_info.image_spacing p.distance_to_crystal = self.user_info.distance_to_crystal @@ -226,13 +232,14 @@ def apply_with_lock(self, actor): # build the data x = np.column_stack((px, py, theta, phi, energy)) + self.debug_nb_hits_before += len(x) # apply the neural network if self.user_info.verbose_batch: print( f"Apply ARF to {energy.shape[0]} hits (device = {self.model_data['current_gpu_mode']})" ) - ax = x[:, 2:5] # two angles and energy + ax = x[:, 2:5] # two angles and energy # FIXME index ? w = self.garf.nn_predict(self.model, self.nn["model_data"], ax) # positions @@ -249,7 +256,7 @@ def apply_with_lock(self, actor): # do nothing if there is no hit in the image if u.shape[0] != 0: - temp = np.zeros(p.image_size, dtype=np.float64) + """temp = np.zeros(p.image_size, dtype=np.float64) temp = self.garf.image_from_coordinates(temp, u, v, w_pred) # add to previous, at the correct slice location # the slice is : current_ene_window + run_id * nb_ene_windows @@ -258,7 +265,14 @@ def apply_with_lock(self, actor): s = p.nb_ene * run_id self.output_image[s : s + p.nb_ene] = ( self.output_image[s : s + p.nb_ene] + temp - ) + )""" + + # alternative + run_id = actor.GetCurrentRunId() + s = p.nb_ene * run_id + img = self.output_image[s : s + p.nb_ene] + self.garf.image_from_coordinates_add(img, u, v, w_pred) + self.debug_nb_hits += u.shape[0] def EndSimulationAction(self): g4.GateARFActor.EndSimulationAction(self) @@ -295,3 +309,7 @@ def EndSimulationAction(self): # write ? if self.user_info.output: itk.imwrite(self.output_image, check_filename_type(self.user_info.output)) + + # FIXME debug + print(f"{self.debug_nb_hits_before=}") + print(f"{self.debug_nb_hits=}") diff --git a/opengate/bin/opengate_tests b/opengate/bin/opengate_tests index 873010bb6..c72e76a6a 100755 --- a/opengate/bin/opengate_tests +++ b/opengate/bin/opengate_tests @@ -56,7 +56,7 @@ def go(test_id, random_tests): except: torch = False - windowsWrongTests = [ + ignored_tests = [ # "test014_engine_2.py", # "test060_PhsSource_ParticleName_direct.py", # "test060_PhsSource_ParticleName_fromPHS_PDGCode.py", @@ -66,6 +66,9 @@ def go(test_id, random_tests): # "test043_garf5_region_MT_subproc.py", # "test043_garf2_region_subproc.py.py", # "test061_user_event_info.py", + "test045_speedup", # this is a binary (still work in progress) + "test066_spect_gaga_garf_2.py", # ignored because reference data (too long) + "test066_spect_gaga_garf_1_reference.py", # ignored because reference data (too long) ] onlyfiles = [ @@ -93,13 +96,11 @@ def go(test_id, random_tests): continue if "_base" in f: continue - if "test045_speedup" in f: - continue if "_helpers" in f: continue if os.name == "nt" and "_mt" in f: continue - if os.name == "nt" and f in windowsWrongTests: + if f in ignored_tests: continue if not torch and f in torch_tests: print(f"Ignoring: {f:<40} (Torch is not available) ") diff --git a/opengate/contrib/spect/genm670.py b/opengate/contrib/spect/genm670.py index df7f2e902..147263788 100644 --- a/opengate/contrib/spect/genm670.py +++ b/opengate/contrib/spect/genm670.py @@ -1,6 +1,7 @@ import pathlib from opengate.exception import fatal from opengate.utility import g4_units +from opengate.managers import Simulation from opengate.geometry.volumes import RepeatParametrisedVolume, HexagonVolume from opengate.geometry.utility import ( @@ -422,7 +423,7 @@ def add_digitizer_energy_windows(sim, crystal_volume_name, channels): return cc -def get_volume_position_in_head(sim, spect_name, vol_name, pos="max"): +def get_volume_position_in_head(sim, spect_name, vol_name, pos="max", axis=2): vol = sim.volume_manager.volumes[f"{spect_name}_{vol_name}"] pMin, pMax = vol.bounding_limits x = pMax @@ -432,12 +433,22 @@ def get_volume_position_in_head(sim, spect_name, vol_name, pos="max"): x = pMin + (pMax - pMin) / 2.0 x = vec_g4_as_np(x) x = translate_point_to_volume(sim, vol, spect_name, x) - return x[2] + return x[axis] + + +def compute_plane_position_and_distance_to_crystal(collimator_type): + sim = Simulation() + spect, crystal = add_ge_nm67_spect_head(sim, "spect", collimator_type, debug=True) + pp = get_volume_position_in_head(sim, "spect", "collimator_psd", "max") + y = get_volume_position_in_head(sim, "spect", "crystal", "center") + dc = pp - y + psd = spect.size[2] / 2.0 - pp + return pp, dc, psd def get_plane_position_and_distance_to_crystal(collimator_type): """ - This has been computed with t043_distances + This has been computed with t043_distances or compute_plane_position_and_distance_to_crystal - first : distance from head center to the PSD (translation for the plane) - second: distance from PSD to center of the crystal - third : distance from the head boundary to the PSD (for spect_radius info) diff --git a/opengate/tests/data b/opengate/tests/data index 19376c4fa..4347b2502 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit 19376c4fa97c41c20f4b96b1e946cff0a1e5bfee +Subproject commit 4347b25021125285a7edbd9ad7d555be9f23e93c diff --git a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py index a7e8609d8..9ee8ba9e4 100644 --- a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py +++ b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py @@ -236,9 +236,7 @@ def gen_cond(n): gsource, 210 * mm, gen_cond ) gsource.generator = gen - gsource.gpu_mode = ( - utility.get_gpu_mode() - ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors + gsource.gpu_mode = utility.get_gpu_mode() def add_gaga_source_vox_condition(sim, p): diff --git a/opengate/tests/utility.py b/opengate/tests/utility.py index 5670c3b4a..1509aa5c4 100644 --- a/opengate/tests/utility.py +++ b/opengate/tests/utility.py @@ -10,10 +10,10 @@ import uproot import sys import matplotlib.pyplot as plt - +from matplotlib.ticker import StrMethodFormatter import gatetools.phsp as phsp +from matplotlib.patches import Circle -# from .helpers_log import colorlog from ..utility import g4_units, check_filename_type from ..exception import fatal, color_error, color_ok from ..image import get_info_from_image, itk_image_view_from_array @@ -252,7 +252,7 @@ def assert_images( sum_tolerance=5, scaleImageValuesFactor=None, ): - # read image and info (size, spacing etc) + # read image and info (size, spacing, etc.) ref_filename1 = check_filename_type(ref_filename1) filename2 = check_filename_type(filename2) img1 = itk.imread(ref_filename1) @@ -1559,3 +1559,244 @@ def get_gpu_mode(): print("Detection of Github actions and MacOS -> Use CPU") return "cpu" return "auto" + + +def np_img_window_level(img, window_width, window_level): + """ + Clip and rescale the grey level values of the image according to window width/level. + Output image is within the range [0, 1] + + Parameters + ---------- + img an image as a numpy array + window_width + window_level + + Returns the clipped and normalized image (range [0, 1]) + ------- + """ + # Apply window/level adjustment to the images + window_min = window_level - window_width / 2 + window_max = window_level + window_width / 2 + clipped_image = np.clip(img, window_min, window_max) + + # Normalize the intensities to the range [0, 1] + normalized_image = (clipped_image - window_min) / (window_max - window_min) + + return normalized_image + + +def np_img_crop(img, crop_center, crop_width): + c = crop_center + w = crop_width + x1 = c[0] - w[0] // 2 + x2 = c[0] + w[0] // 2 + y1 = c[1] - w[1] // 2 + y2 = c[1] + w[1] // 2 + img = img[:, y1:y2, x1:x2] + return img, (x1, x2, y1, y2) + + +def np_plot_slice( + ax, + img, + num_slice, + window_width, + window_level, + crop_center, + crop_width, + spacing=(1, 1), +): + # crop and grey level + img = np_img_window_level(img, window_width, window_level) + img, crop_coord = np_img_crop(img, crop_center, crop_width) + + # slice + slice = img[num_slice, :, :] + im = ax.imshow(slice, cmap="gray") + + # prepare ticks + nticks = 6 + x_step = int(np.around((crop_coord[1] - crop_coord[0]) / nticks)) + x_ticks = np.char.mod( + "%.0f", + np.around( + np.arange(crop_coord[0], crop_coord[1], x_step) * spacing[0], decimals=1 + ), + ) + y_step = int(np.around((crop_coord[3] - crop_coord[2]) / nticks)) + y_ticks = np.char.mod( + "%.0f", + np.around( + np.arange(crop_coord[2], crop_coord[3], y_step) * spacing[1], decimals=1 + ), + ) + + # ticks + ax.set_xticks(np.arange(0, crop_width[0], x_step), x_ticks) + ax.set_yticks(np.arange(0, crop_width[1], y_step), y_ticks) + ax.set_xlabel("X (mm)") + ax.set_ylabel("Y (mm)") + return im + + +def np_plot_slice_h_line(ax, hline, crop_center, crop_width): + x = np.arange(0, 100) + c = int(hline - (crop_center[1] - crop_width[1] / 2)) + y = [c] * len(x) + ax.plot(x, y, color="r") + + +def np_plot_slice_v_line(ax, vline, crop_center, crop_width): + x = np.arange(0, 100) + c = int(vline - (crop_center[0] - crop_width[0] / 2)) + y = [c] * len(x) + ax.plot(y, x, color="r") + + +def add_colorbar(imshow, window_level, window_width): + cbar = plt.colorbar( + imshow, orientation="vertical", format=StrMethodFormatter("{x:.1f}") + ) + # window_min = window_level - window_width / 2 + window_max = window_level + window_width / 2 + # Number of ticks you want on the color bar + num_ticks = 10 + tick_values = np.linspace(0, window_max, num_ticks) + cbar.set_ticks(tick_values) + + +def np_plot_integrated_profile( + ax, img, axis, num_slice, crop_center, crop_width, label, spacing +): + img, crop_coord = np_img_crop(img, crop_center, crop_width) + img = img[num_slice, :, :] + profile = np.mean(img, axis=axis) + values = np.arange(0, len(profile)) * spacing + crop_coord[axis * 2] * spacing + ax.plot(values, profile, label=label) + + +def np_plot_profile_X(ax, img, hline, num_slice, crop_center, crop_width, label, width): + c = int(hline - (crop_center[1] - crop_width[1] / 2)) + img, _ = np_img_crop(img, crop_center, crop_width) + if width == 0: + img = img[num_slice, c : c + 1, :] + else: + img = img[num_slice, c - width : c + width, :] + y = np.mean(img, axis=0) + x = np.arange(0, len(y)) + ax.plot(x, y, label=label) + + +def np_plot_profile_Y(ax, img, vline, num_slice, crop_center, crop_width, label, width): + c = int(vline - (crop_center[0] - crop_width[0] / 2)) + img, _ = np_img_crop(img, crop_center, crop_width) + if width == 0: + img = img[num_slice, :, c : c + 1] + else: + img = img[num_slice, :, c - width : c + width] + x = np.mean(img, axis=1) + y = np.arange(0, len(x)) + ax.plot(y, x, label=label) + + +def np_get_circle_mean_value(img, center, radius): + y, x = np.ogrid[: img.shape[0], : img.shape[1]] + distance_squared = (x - center[0]) ** 2 + (y - center[1]) ** 2 + mask = distance_squared <= radius**2 + pixels_within_circle = img[mask] + mean_value = np.mean(pixels_within_circle) + return mean_value + + +def add_circle(ax, img, crop_center, crop_width, center, radius): + _, crop = np_img_crop(img, crop_center, crop_width) + circle = Circle( + (center[0] - crop[0], center[1] - crop[2]), + radius, + linewidth=2, + edgecolor="r", + facecolor="none", + ) + ax.add_patch(circle) + + +def add_border(ax, border_color, border_width): + # Set the spines color and width + for spine in ax.spines.values(): + spine.set_edgecolor(border_color) + spine.set_linewidth(border_width) + + +def plot_compare_profile(ref_names, test_names, options): + # options + scaling = options.scaling + n_slice = options.n_slice + ww = options.window_width + wl = options.window_level + c = options.crop_center + w = options.crop_width + hline = options.hline + vline = options.vline + wi = options.width + lab_ref = options.lab_ref + lab_test = options.lab_test + title = options.title + + # read as np array + img_ref = [] + img_test = [] + for ref_name, test_name in zip(ref_names, test_names): + iref = itk.imread(ref_name) + spacing = (iref.GetSpacing()[1], iref.GetSpacing()[2]) + iref = itk.array_from_image(iref) + itest = itk.imread(test_name) + itest = itk.array_from_image(itest) * scaling + img_ref.append(iref) + img_test.append(itest) + + # plot + n = len(img_ref) + nrow = 2 + ncol = 2 * n + _, ax = plt.subplots(nrow, ncol, figsize=(ncol * 6, 10)) + for i in range(n): + np_plot_slice(ax[0][i * n], img_ref[i], n_slice, ww, wl, c, w, spacing) + last = np_plot_slice( + ax[0][i * n + 1], img_test[i], n_slice, ww, wl, c, w, spacing + ) + np_plot_slice_h_line(ax[0][i * n], hline, c, w) + np_plot_slice_h_line(ax[0][i * n + 1], hline, c, w) + np_plot_slice_v_line(ax[0][i * n], vline, c, w) + np_plot_slice_v_line(ax[0][i * n + 1], vline, c, w) + + # Add colorbar to the figure + add_colorbar(last, wl, ww) + + # profiles + lref = f"{lab_ref} (horizontal)" + ltest = f"{lab_test} (horizontal)" + for i in range(len(img_ref)): + np_plot_profile_X( + ax[1][i * n], img_ref[i], hline, n_slice, c, w, lref, width=wi + ) + np_plot_profile_X( + ax[1][i * n], img_test[i], hline, n_slice, c, w, ltest, width=wi + ) + ax[1][i * n].legend() + + lref = f"{lab_ref} (vertical)" + ltest = f"{lab_test} (vertical)" + for i in range(len(img_ref)): + np_plot_profile_Y( + ax[1][i * n + 1], img_ref[i], vline, n_slice, c, w, lref, width=wi + ) + np_plot_profile_Y( + ax[1][i * n + 1], img_test[i], vline, n_slice, c, w, ltest, width=wi + ) + ax[1][i * n + 1].legend() + + plt.suptitle(title, fontweight="bold", fontsize=12, color="red") + # Adjust spacing between subplots if necessary + plt.tight_layout() + return plt From 21e0e6d3f3fb30aea0490dfa7e0b6c7448b91fbe Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 14 Nov 2023 11:46:24 +0100 Subject: [PATCH 06/37] adapt to new garf gaga --- opengate/actors/arfactors.py | 16 +++-------- opengate/sources/gansources.py | 27 +++++-------------- opengate/tests/src/test045_speedup.py | 5 +++- .../test066_spect_gaga_garf_1_reference.py | 18 +++++++++++++ .../tests/src/test066_spect_gaga_garf_2.py | 3 +++ .../test066_spect_gaga_garf_3_standalone.py | 15 +++++++++-- .../src/test066_spect_gaga_garf_4_analyse1.py | 17 +++++++----- .../src/test066_spect_gaga_garf_5_analyse2.py | 19 +++++++------ opengate/tests/utility.py | 4 +-- 9 files changed, 70 insertions(+), 54 deletions(-) diff --git a/opengate/actors/arfactors.py b/opengate/actors/arfactors.py index 87d9391ec..edd2f1e85 100644 --- a/opengate/actors/arfactors.py +++ b/opengate/actors/arfactors.py @@ -252,22 +252,12 @@ def apply_with_lock(self, actor): coord = np.around(coord).astype(int) v = coord[:, 0] u = coord[:, 1] - u, v, w_pred = self.garf.remove_out_of_image_boundaries(u, v, w, p.image_size) + u, v, w_pred = self.garf.remove_out_of_image_boundaries2( + u, v, w, self.user_info.image_size + ) # do nothing if there is no hit in the image if u.shape[0] != 0: - """temp = np.zeros(p.image_size, dtype=np.float64) - temp = self.garf.image_from_coordinates(temp, u, v, w_pred) - # add to previous, at the correct slice location - # the slice is : current_ene_window + run_id * nb_ene_windows - run_id = actor.GetCurrentRunId() - # self.simulation_engine_wr().g4_RunManager.GetCurrentRun().GetRunID() - s = p.nb_ene * run_id - self.output_image[s : s + p.nb_ene] = ( - self.output_image[s : s + p.nb_ene] + temp - )""" - - # alternative run_id = actor.GetCurrentRunId() s = p.nb_ene * run_id img = self.output_image[s : s + p.nb_ene] diff --git a/opengate/sources/gansources.py b/opengate/sources/gansources.py index 3e85f9bda..c4770081a 100644 --- a/opengate/sources/gansources.py +++ b/opengate/sources/gansources.py @@ -427,7 +427,7 @@ def read_gan_and_keys(self): self.gan_info = Box() g = self.gan_info g.params, g.G, g.D, g.optim = self.gaga.load( - self.user_info.pth_filename, self.gpu_mode, verbose=False + self.user_info.pth_filename, self.gpu_mode ) """ @@ -576,10 +576,9 @@ def generator(self, source): print(f"Generate {n} particles from GAN ", end="") # generate samples (this is the most time-consuming part) - fake = self.gaga.generate_samples2( + fake = self.gaga.generate_samples_no_cond( g.params, g.G, - g.D, n=n, batch_size=n, normalize=False, @@ -747,10 +746,9 @@ def generator(self, source): print(f"Generate {n} particles from GAN ", end="") # generate samples (this is the most time-consuming part) - fake = self.gaga.generate_samples2( + fake = self.gaga.generate_samples_no_cond( g.params, g.G, - g.D, n=n, batch_size=n, normalize=False, @@ -896,16 +894,11 @@ def generator(self, source): # needed by test 047 fake = cond else: - fake = self.gaga.generate_samples2( + fake = self.gaga.generate_samples3( g.params, g.G, - g.D, n=n, - batch_size=n, - normalize=False, - to_numpy=True, cond=cond, - silence=True, ) # consider the names of the output keys position/direction/energy/time/weight @@ -975,16 +968,8 @@ def generator(self, source): cond = self.generate_condition(n) # generate samples (this is the most time-consuming part) - fake = self.gaga.generate_samples2( - g.params, - g.G, - g.D, - n=n, - batch_size=n, - normalize=False, - to_numpy=False, # next step (from_tlor_to_pairs) is in torch, not numpy - cond=cond, - silence=True, + fake = self.gaga.generate_samples3( + g.params, g.G, to_numpy=False, n=n, cond=cond ) # parametrisation diff --git a/opengate/tests/src/test045_speedup.py b/opengate/tests/src/test045_speedup.py index 3778eeb10..7e4f1c2e7 100755 --- a/opengate/tests/src/test045_speedup.py +++ b/opengate/tests/src/test045_speedup.py @@ -22,7 +22,10 @@ @click.option("--threads", "-t", default=1, help="Nb of threads") @click.option("--visu", "-v", default=False, is_flag=True, help="visu for debug") @click.option( - "--output_folder", "-o", default=".", help="output folder (AUTO for the tests)" + "--output_folder", + "-o", + default=".", + help="output folder (AUTO for the tests)", ) @click.option("--seed", default="auto", help="random engine seed") def go( diff --git a/opengate/tests/src/test066_spect_gaga_garf_1_reference.py b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py index 7f2452169..2afe21db6 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_1_reference.py +++ b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py @@ -9,6 +9,24 @@ paths = utility.get_default_test_paths(__file__, None, "test066") output_path = paths.output + """ + Test066 + Spect image simulation of IEC phantom with Tc99m source, two heads (two projections). + With analog methods and gaga+garf method (within and outside Gate). + The activity concentration is the same in all spheres, 100 kBq, so 4.78 MBq in total. + No background activity. + The total events for 30 sec acquisition (no half life) is around 1.3e8 gammas. + + Test on linux + - test066_1 : reference + - test066_2 : gaga + garf, within gate (around ~x4 faster than reference, and need less particles) + - test066_3 : gaga + garf, outside gate, only one thread, ~50% faster than within gate + + test066_2 uses 2x less particles than test066_1 (because ARF) + test066_3 uses 2x less particles than test066_2 (to save time) + + """ + # create the simulation sim = gate.Simulation() simu_name = "test066_1_reference" diff --git a/opengate/tests/src/test066_spect_gaga_garf_2.py b/opengate/tests/src/test066_spect_gaga_garf_2.py index abf22d43e..af18d8162 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_2.py +++ b/opengate/tests/src/test066_spect_gaga_garf_2.py @@ -12,10 +12,12 @@ # create the simulation sim = gate.Simulation() simu_name = "test066_2" + simu_name = "test066_2_fake" # FIXME # options ui = sim.user_info ui.number_of_threads = 5 + ui.number_of_threads = 1 # FIXME # ui.visu = True ui.visu_type = "vrml" ui.random_seed = "auto" @@ -34,6 +36,7 @@ # Estimated gammas 127,008,708 gammas (weights = 0.8850) # so, because we use ARF, about 1/2 particles needed total_activity = 127008708 / 30 * Bq / ui.number_of_threads / 2 + total_activity = 127008708 / 30 * Bq / ui.number_of_threads / 40 # FIXME print(f"Total activity: {total_activity/Bq:.0f} Bq") if ui.visu: total_activity = 10 * Bq diff --git a/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py b/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py index 44e540f92..17eb29e42 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py +++ b/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py @@ -8,6 +8,7 @@ from scipy.spatial.transform import Rotation import opengate as gate from opengate.contrib.spect import genm670 +import time if __name__ == "__main__": paths = utility.get_default_test_paths(__file__, None, "test066") @@ -37,6 +38,7 @@ "batch_size": 1e5, "gpu_mode": "auto", "backward_distance": 50 * mm, + "verbose": 0, } ) garf_user_info = Box( @@ -48,6 +50,8 @@ "distance_to_crystal": crystal_dist, "batch_size": 1e5, "gpu_mode": "auto", + "verbose": 0, + "hit_slice": False, } ) @@ -60,15 +64,22 @@ r = Rotation.from_euler("y", 0, degrees=True) garf_user_info.plane_rotation = r - # GO - n = 127008708 / 20 + # define the angles + n = 127008708 / 4 angle_rotations = [ Rotation.from_euler("y", 0, degrees=True), Rotation.from_euler("y", 180, degrees=True), ] + + # GO + t1 = time.time() images = gaga.gaga_garf_generate_spect( gaga_user_info, garf_user_info, n, angle_rotations ) + t2 = time.time() + t = t2 - t1 + print(f"Computation time is {t:.2f} seconds") + print(f"Computation PPS is {n/t:.0f} ") # save image i = 0 diff --git a/opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py b/opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py index 813c0a2f2..118e9ecc0 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py +++ b/opengate/tests/src/test066_spect_gaga_garf_4_analyse1.py @@ -48,13 +48,16 @@ print() is_ok = True for ref_name, test_name in zip(ref_names, test_names): - is_ok = utility.assert_images( - ref_name, - test_name, - axis="x", - scaleImageValuesFactor=options.scaling, - sum_tolerance=8, - tolerance=110, + is_ok = ( + utility.assert_images( + ref_name, + test_name, + axis="x", + scaleImageValuesFactor=options.scaling, + sum_tolerance=12, + tolerance=110, + ) + and is_ok ) utility.test_ok(is_ok) diff --git a/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py b/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py index 0c77182ae..8566290bf 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py +++ b/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py @@ -25,7 +25,7 @@ ] options = Box() - options.scaling = 10 + options.scaling = 2 options.n_slice = 1 options.window_width = 33 options.window_level = 17 @@ -48,13 +48,16 @@ print() is_ok = True for ref_name, test_name in zip(ref_names, test_names): - is_ok = utility.assert_images( - ref_name, - test_name, - axis="x", - scaleImageValuesFactor=options.scaling, - sum_tolerance=15, - tolerance=105, + is_ok = ( + utility.assert_images( + ref_name, + test_name, + axis="x", + scaleImageValuesFactor=options.scaling, + sum_tolerance=15, + tolerance=120, + ) + and is_ok ) utility.test_ok(is_ok) diff --git a/opengate/tests/utility.py b/opengate/tests/utility.py index 1509aa5c4..3e005e3d7 100644 --- a/opengate/tests/utility.py +++ b/opengate/tests/utility.py @@ -1749,9 +1749,9 @@ def plot_compare_profile(ref_names, test_names, options): for ref_name, test_name in zip(ref_names, test_names): iref = itk.imread(ref_name) spacing = (iref.GetSpacing()[1], iref.GetSpacing()[2]) - iref = itk.array_from_image(iref) + iref = itk.array_view_from_image(iref) itest = itk.imread(test_name) - itest = itk.array_from_image(itest) * scaling + itest = itk.array_view_from_image(itest) * scaling img_ref.append(iref) img_test.append(itest) From 31180f0eec38bbcb304147207a525b98675e8fa8 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 14 Nov 2023 11:48:25 +0100 Subject: [PATCH 07/37] ignore t66 --- opengate/bin/opengate_tests | 1 + 1 file changed, 1 insertion(+) diff --git a/opengate/bin/opengate_tests b/opengate/bin/opengate_tests index 7fdf5d58e..47e706e91 100755 --- a/opengate/bin/opengate_tests +++ b/opengate/bin/opengate_tests @@ -69,6 +69,7 @@ def go(test_id, random_tests): "test045_speedup", # this is a binary (still work in progress) "test066_spect_gaga_garf_2.py", # ignored because reference data (too long) "test066_spect_gaga_garf_1_reference.py", # ignored because reference data (too long) + "test066_spect_gaga_garf_3_standalone.py", # ignored becausetoo long ] onlyfiles = [ From 1155be49cca5e5a3fcf690804cfef83537a80fba Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 17 Nov 2023 09:21:32 +0100 Subject: [PATCH 08/37] make AA work for parallel world --- .../opengate_lib/GateHelpersGeometry.cpp | 6 +- opengate/bin/opengate_tests | 18 +-- opengate/contrib/spect/genm670.py | 142 +++++++++++++++--- .../test066_spect_gaga_garf_1_reference.py | 4 +- .../tests/src/test066_spect_gaga_garf_2.py | 9 +- .../test066_spect_gaga_garf_3_standalone.py | 14 +- .../src/test066_spect_gaga_garf_helpers.py | 91 +++++------ 7 files changed, 195 insertions(+), 89 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateHelpersGeometry.cpp b/core/opengate_core/opengate_lib/GateHelpersGeometry.cpp index 875b8486d..4febe8557 100644 --- a/core/opengate_core/opengate_lib/GateHelpersGeometry.cpp +++ b/core/opengate_core/opengate_lib/GateHelpersGeometry.cpp @@ -24,7 +24,11 @@ void ComputeTransformationFromVolumeToWorld(const std::string &phys_volume_name, auto rot = phys->GetObjectRotationValue(); rotation = rot * rotation; translation = rot * translation + tr; - name = phys->GetMotherLogical()->GetName(); + // Warning, the world can be a parallel world + if (phys->GetMotherLogical() == nullptr) + name = "world"; + else + name = phys->GetMotherLogical()->GetName(); } } diff --git a/opengate/bin/opengate_tests b/opengate/bin/opengate_tests index 47e706e91..5147602cc 100755 --- a/opengate/bin/opengate_tests +++ b/opengate/bin/opengate_tests @@ -33,7 +33,7 @@ def go(test_id, random_tests): pathlib.Path(opengate.tests.__file__).resolve().parent, "../tests/src" ) - print("Look for tests in: " + mypath) + print("Looking for tests in: " + mypath) if not check_tests_data_folder(): return False @@ -57,19 +57,11 @@ def go(test_id, random_tests): torch = False ignored_tests = [ - # "test014_engine_2.py", - # "test060_PhsSource_ParticleName_direct.py", - # "test060_PhsSource_ParticleName_fromPHS_PDGCode.py", - # "test060_PhsSource_ParticleName_fromPHS_ParticleName.py", - # "test060_PhsSource_rotation.py", - # "test060_PhsSource_translation.py", - # "test043_garf5_region_MT_subproc.py", - # "test043_garf2_region_subproc.py.py", - # "test061_user_event_info.py", "test045_speedup", # this is a binary (still work in progress) - "test066_spect_gaga_garf_2.py", # ignored because reference data (too long) + "test066_spect_gaga_garf_0_orientation.py", # ignored because visu only "test066_spect_gaga_garf_1_reference.py", # ignored because reference data (too long) - "test066_spect_gaga_garf_3_standalone.py", # ignored becausetoo long + "test066_spect_gaga_garf_2.py", # ignored because reference data (too long, GPU) + "test066_spect_gaga_garf_3_standalone.py", # ignored because too long (GPU) ] onlyfiles = [ @@ -79,7 +71,7 @@ def go(test_id, random_tests): files = [] for f in onlyfiles: if "wip" in f: - print(f"Ignoring: {f:<40} (Work In Progress) ") + print(f"Ignoring: {f:<40} ") continue if "visu" in f: continue diff --git a/opengate/contrib/spect/genm670.py b/opengate/contrib/spect/genm670.py index 147263788..f9af73b0c 100644 --- a/opengate/contrib/spect/genm670.py +++ b/opengate/contrib/spect/genm670.py @@ -2,26 +2,14 @@ from opengate.exception import fatal from opengate.utility import g4_units from opengate.managers import Simulation - from opengate.geometry.volumes import RepeatParametrisedVolume, HexagonVolume from opengate.geometry.utility import ( translate_point_to_volume, get_transform_orbiting, vec_g4_as_np, ) - -# unit -cm = g4_units.cm -mm = g4_units.mm -deg = g4_units.deg - -# colors -red = [1, 0.7, 0.7, 0.8] -blue = [0.5, 0.5, 1, 0.8] -gray = [0.5, 0.5, 0.5, 1] -white = [1, 1, 1, 1] -yellow = [1, 1, 0, 1] -green = [0, 1, 0, 1] +from scipy.spatial.transform import Rotation +from box import Box def get_collimator(rad): @@ -33,6 +21,8 @@ def get_collimator(rad): def add_ge_nm67_fake_spect_head(sim, name="spect"): + white = [1, 1, 1, 1] + cm = g4_units.cm spect_length = 19 * cm head = sim.add_volume("Box", name) head.material = "G4_AIR" @@ -78,7 +68,7 @@ def add_ge_nm67_spect_head(sim, name="spect", collimator_type="lehr", debug=Fals # spect collimator if collimator_type: - colli = add_ge_nm670_spect_collimator(sim, name, head, collimator_type, debug) + add_ge_nm670_spect_collimator(sim, name, head, collimator_type, debug) return head, crystal @@ -93,6 +83,15 @@ def distance_to_center_of_crystal(sim, name="spect"): def add_ge_nm670_spect_box(sim, name, collimator_type): + cm = g4_units.cm + + # colors + blue = [0.5, 0.5, 1, 0.8] + gray = [0.5, 0.5, 0.5, 1] + white = [1, 1, 1, 1] + yellow = [1, 1, 0, 1] + green = [0, 1, 0, 1] + # the total length spect_length = 19 * cm @@ -149,6 +148,8 @@ def add_ge_nm670_spect_box(sim, name, collimator_type): def add_ge_nm670_spect_crystal(sim, name, lead_cover): + cm = g4_units.cm + yellow = [1, 1, 0, 1] # mono-bloc crystal thickness 3/8 of inch = 0.9525 cm # (if 5/8 inch = 1.5875 ; but probably need to translate elements) crystal = sim.add_volume("Box", f"{name}_crystal") @@ -165,6 +166,12 @@ def add_ge_nm670_spect_collimator(sim, name, head, collimator_type, debug): Start with default lehr collimator description, then change some parameters for the other types """ + cm = g4_units.cm + + # colors + red = [1, 0.7, 0.7, 0.8] + blue = [0.5, 0.5, 1, 0.8] + green = [0, 1, 0, 1] # mono-bloc crystal thickness 3/8 of inch colli_trd = sim.add_volume("Trd", f"{name}_collimator_trd") @@ -256,6 +263,8 @@ def add_ge_nm670_spect_collimator(sim, name, head, collimator_type, debug): def hegp_collimator_repeater(sim, name, core, debug): + cm = g4_units.cm + mm = g4_units.mm # one single hole hole = sim.add_volume("Hexagon", f"{name}_collimator_hole") hole.height = 6.6 * cm @@ -277,6 +286,8 @@ def hegp_collimator_repeater(sim, name, core, debug): def megp_collimator_repeater(sim, name, core, debug): + cm = g4_units.cm + mm = g4_units.mm # one single hole hole = sim.add_volume("Hexagon", f"{name}_collimator_hole") hole.height = 5.8 * cm @@ -300,6 +311,8 @@ def megp_collimator_repeater(sim, name, core, debug): def lehr_collimator_repeater(sim, name, core, debug): + cm = g4_units.cm + mm = g4_units.mm # one single hole hole = HexagonVolume(name=f"{name}_collimator_hole") hole.height = 3.5 * cm @@ -423,6 +436,7 @@ def add_digitizer_energy_windows(sim, crystal_volume_name, channels): return cc +# FIXME : put this elsewhere def get_volume_position_in_head(sim, spect_name, vol_name, pos="max", axis=2): vol = sim.volume_manager.volumes[f"{spect_name}_{vol_name}"] pMin, pMax = vol.bounding_limits @@ -439,11 +453,11 @@ def get_volume_position_in_head(sim, spect_name, vol_name, pos="max", axis=2): def compute_plane_position_and_distance_to_crystal(collimator_type): sim = Simulation() spect, crystal = add_ge_nm67_spect_head(sim, "spect", collimator_type, debug=True) - pp = get_volume_position_in_head(sim, "spect", "collimator_psd", "max") + pos = get_volume_position_in_head(sim, "spect", "collimator_psd", "max") y = get_volume_position_in_head(sim, "spect", "crystal", "center") - dc = pp - y - psd = spect.size[2] / 2.0 - pp - return pp, dc, psd + crystal_distance = pos - y + psd = spect.size[2] / 2.0 - pos + return pos, crystal_distance, psd def get_plane_position_and_distance_to_crystal(collimator_type): @@ -465,3 +479,93 @@ def get_plane_position_and_distance_to_crystal(collimator_type): fatal( f'Unknown collimator type "{collimator_type}", please use lehr or megp or hegp' ) + + +def add_fake_table(sim, name="table"): + """ + Add a patient table (fake) + """ + + # unit + mm = g4_units.mm + cm = g4_units.cm + cm3 = g4_units.cm3 + deg = g4_units.deg + gcm3 = g4_units.g / cm3 + + # colors + red = [1, 0.7, 0.7, 0.8] + white = [1, 1, 1, 1] + + sim.add_material_weights(f"CarbonFiber", ["C"], [1], 1.78 * gcm3) + + # main bed + table = sim.add_volume("Tubs", f"{name}_table") + table.mother = "world" + table.rmax = 439 * mm + table.rmin = 406 * mm + table.dz = 200 * cm / 2.0 + table.sphi = 0 * deg + table.dphi = 70 * deg + table.translation = [0, 25 * cm, 0] + table.rotation = Rotation.from_euler("z", -125, degrees=True).as_matrix() + table.material = "CarbonFiber" + table.color = white + + # interior of the table + tablein = sim.add_volume("Tubs", f"{name}_tablein") + tablein.mother = table.name + tablein.rmax = 436.5 * mm + tablein.rmin = 408.5 * mm + tablein.dz = 200 * cm / 2.0 + tablein.sphi = 0 * deg + tablein.dphi = 69 * deg + tablein.translation = [0, 0, 0] + tablein.rotation = Rotation.from_euler("z", 0.5, degrees=True).as_matrix() + tablein.material = "G4_AIR" + tablein.color = red + + return table + + +def set_head_orientation(head, collimator_type, radius, gantry_angle=0): + # pos is the distance from entrance detection plane and head boundary + pos, _, _ = compute_plane_position_and_distance_to_crystal(collimator_type) + distance = radius + pos + # rotation X180 is to set the detector head-foot + # rotation Z90 is the gantry angle + r1 = Rotation.from_euler("x", 90, degrees=True) + r2 = Rotation.from_euler("z", gantry_angle, degrees=True) + r = r2 * r1 + head.translation = r.apply([0, 0, -distance]) + head.rotation = r.as_matrix() + return r + + +def add_detection_plane_for_arf( + sim, plane_size, colli_type, radius, gantry_angle=0, det_name=None +): + if det_name is None: + det_name = "arf_plane" + + # rotation ? no + r = Rotation.from_euler("yx", (0, 0), degrees=True) + + # plane + nm = g4_units.nm + detector_plane = sim.add_volume("Box", det_name) + detector_plane.material = "G4_Galactic" + detector_plane.color = [1, 0, 0, 1] + detector_plane.size = [plane_size[0], plane_size[1], 1 * nm] + + # (fake) initial head rotation + head = Box() + head.translation = None + head.rotation = None + ri = set_head_orientation(head, colli_type, radius, gantry_angle) + + # orientation + detector_plane.rotation = (ri * r).as_matrix() + detector_plane.translation = ri.apply([0, 0, -radius]) + + return detector_plane diff --git a/opengate/tests/src/test066_spect_gaga_garf_1_reference.py b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py index 2afe21db6..836e70e7d 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_1_reference.py +++ b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py @@ -37,7 +37,6 @@ # ui.visu = True ui.visu_type = "vrml" ui.random_seed = "auto" - ui.check_volumes_overlap = False # units mm = gate.g4_units.mm @@ -52,8 +51,9 @@ proj2.output = f"{output_path}/{simu_name}_1.mhd" # add IEC phantom - gate_iec.add_iec_phantom(sim, name="iec") + iec = gate_iec.add_iec_phantom(sim, name="iec") sim.set_production_cut("iec", "all", 1 * mm) + iec.rotation = Rotation.from_euler("x", 90, degrees=True).as_matrix() # sources IEC ac = 1e5 * BqmL # 1e5 = about 10 min with 10 threads linux diff --git a/opengate/tests/src/test066_spect_gaga_garf_2.py b/opengate/tests/src/test066_spect_gaga_garf_2.py index af18d8162..f11464a5d 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_2.py +++ b/opengate/tests/src/test066_spect_gaga_garf_2.py @@ -12,16 +12,13 @@ # create the simulation sim = gate.Simulation() simu_name = "test066_2" - simu_name = "test066_2_fake" # FIXME # options ui = sim.user_info ui.number_of_threads = 5 - ui.number_of_threads = 1 # FIXME # ui.visu = True ui.visu_type = "vrml" ui.random_seed = "auto" - # ui.check_volumes_overlap = False # units mm = gate.g4_units.mm @@ -31,15 +28,15 @@ BqmL = Bq / cm3 # activity - # FIXME to compute exact number + # FIXME compute exact number w, e = get_rad_gamma_energy_spectrum("Tc99m") # Estimated gammas 127,008,708 gammas (weights = 0.8850) # so, because we use ARF, about 1/2 particles needed total_activity = 127008708 / 30 * Bq / ui.number_of_threads / 2 - total_activity = 127008708 / 30 * Bq / ui.number_of_threads / 40 # FIXME + # total_activity = 127008708 / 30 * Bq / ui.number_of_threads / 40 # FIXME print(f"Total activity: {total_activity/Bq:.0f} Bq") if ui.visu: - total_activity = 10 * Bq + total_activity = 1 * Bq # source # voxelize_iec_phantom -o ../data/iec_2mm.mhd -s 2 --output_source data/iec_source_same_spheres_2mm.mhd diff --git a/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py b/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py index 17eb29e42..3ac5d4200 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py +++ b/opengate/tests/src/test066_spect_gaga_garf_3_standalone.py @@ -46,7 +46,7 @@ "pth_filename": garf_pth_filename, "image_size": [128, 128], "image_spacing": [3 * mm, 3 * mm], - "plane_distance": -head_radius - pos, + "plane_distance": head_radius, "distance_to_crystal": crystal_dist, "batch_size": 1e5, "gpu_mode": "auto", @@ -60,18 +60,22 @@ # initialize gaga and garf (read the NN) gaga.gaga_garf_generate_spect_initialize(gaga_user_info, garf_user_info) - # Initial rotation angle of the head (identity here) - r = Rotation.from_euler("y", 0, degrees=True) + # Initial rotation of the iec -> X90 inverted + r_iec = Rotation.from_euler("x", -90, degrees=True) + + # Initial rotation angle of the head + r = Rotation.from_euler("x", 90, degrees=True) + r = r * r_iec garf_user_info.plane_rotation = r # define the angles - n = 127008708 / 4 angle_rotations = [ Rotation.from_euler("y", 0, degrees=True), - Rotation.from_euler("y", 180, degrees=True), + Rotation.from_euler("y", 180, degrees=True), # FIXME why Y ????? ] # GO + n = 127008708 / 4 t1 = time.time() images = gaga.gaga_garf_generate_spect( gaga_user_info, garf_user_info, n, angle_rotations diff --git a/opengate/tests/src/test066_spect_gaga_garf_helpers.py b/opengate/tests/src/test066_spect_gaga_garf_helpers.py index ffe83a3b3..dd2e7f557 100644 --- a/opengate/tests/src/test066_spect_gaga_garf_helpers.py +++ b/opengate/tests/src/test066_spect_gaga_garf_helpers.py @@ -12,7 +12,7 @@ def create_world(sim): # world size m = gate.g4_units.m - sim.world.size = [1 * m, 1 * m, 1 * m] + sim.world.size = [2 * m, 2 * m, 2 * m] sim.world.material = "G4_AIR" @@ -38,11 +38,10 @@ def create_simu_with_genm670(sim, collimator_type="lehr", debug=False): sim, "spect2", collimator_type=collimator_type, debug=debug ) - # default rotation to be in front of the IEC (not a realistic orientation) - head1.translation = [0, 0, -280 * mm] - head2.translation = [0, 0, 280 * mm] - head1.rotation = Rotation.from_euler("y", 0, degrees=True).as_matrix() - head2.rotation = Rotation.from_euler("y", 180, degrees=True).as_matrix() + # default rotation to be in front of the IEC + radius = 280 * mm + genm670.set_head_orientation(head1, collimator_type, radius, 0) + genm670.set_head_orientation(head2, collimator_type, radius, 180) # digitizer keV = gate.g4_units.keV @@ -73,12 +72,35 @@ def create_simu_with_gaga( create_world(sim) # add gaga source - add_gaga_source(sim, total_activity, activity_source, gaga_pth_filename) + gsource = add_gaga_source(sim, total_activity, activity_source, gaga_pth_filename) + gsource.position.rotation = Rotation.from_euler("x", 90, degrees=True).as_matrix() # add arf plane mm = gate.g4_units.mm - det1, det2 = add_garf_detector_planes(sim, head_radius=280 * mm) - arf1, arf2 = add_garf_detectors(sim, garf_pth_filename, det1, det2) + colli_type = "lehr" + radius = 280 * mm + size = [128, 128] + spacing = [3 * mm, 3 * mm] + plane_size = [size[0] * spacing[0], size[1] * spacing[1]] + _, crystal_dist, _ = genm670.compute_plane_position_and_distance_to_crystal( + colli_type + ) + det1 = genm670.add_detection_plane_for_arf( + sim, plane_size, colli_type, radius, 0, "det1" + ) + det2 = genm670.add_detection_plane_for_arf( + sim, plane_size, colli_type, radius, 180, "det2" + ) + + print(f"{crystal_dist=}") + + # add actors + arf1 = add_arf_actor( + sim, det1, size, spacing, crystal_dist, "arf1", garf_pth_filename + ) + arf2 = add_arf_actor( + sim, det2, size, spacing, crystal_dist, "arf2", garf_pth_filename + ) # stats sim.add_actor("SimulationStatisticsActor", "stats") @@ -120,15 +142,12 @@ def add_gaga_source(sim, total_activity, activity_source, gaga_pth_filename): gsource, cond_generator.generate_condition ) gsource.generator = gene + return gsource -def add_garf_detector_planes(sim, head_radius): +def add_garf_detector_planes_OLD(sim, plane_size, head_radius): cm = gate.g4_units.cm nm = gate.g4_units.nm - # GE colli size - colli_size = [54.6 * cm, 40.6 * cm] # FIXME function - # colli_size = [100 * cm, 100 * cm] # FIXME function - colli_size = [128 * 3, 128 * 3] # FIXME function pos, crystal_dist, psd = genm670.compute_plane_position_and_distance_to_crystal( "lehr" ) @@ -139,7 +158,7 @@ def add_garf_detector_planes(sim, head_radius): detector_plane1 = sim.add_volume("Box", "det_plane1") detector_plane1.material = "G4_Galactic" detector_plane1.color = [1, 0, 0, 1] - detector_plane1.size = [colli_size[0], colli_size[1], 1 * nm] + detector_plane1.size = [plane_size[0], plane_size[1], 1 * nm] r = Rotation.from_euler("x", 0, degrees=True) detector_plane1.rotation = r.as_matrix() detector_plane1.translation = [0, 0, -head_radius - pos] @@ -148,7 +167,7 @@ def add_garf_detector_planes(sim, head_radius): detector_plane2 = sim.add_volume("Box", "det_plane2") detector_plane2.material = "G4_Galactic" detector_plane2.color = [1, 0, 0, 1] - detector_plane2.size = [colli_size[0], colli_size[1], 1 * nm] + detector_plane2.size = [plane_size[0], plane_size[1], 1 * nm] r = Rotation.from_euler("y", 180, degrees=True) detector_plane2.rotation = r.as_matrix() detector_plane2.translation = [0, 0, -head_radius - pos] @@ -159,33 +178,19 @@ def add_garf_detector_planes(sim, head_radius): return detector_plane1, detector_plane2 -def add_garf_detectors(sim, pth_filename, det1, det2): - mm = gate.g4_units.mm - detectors = [det1, det2] - - pos, crystal_dist, psd = genm670.compute_plane_position_and_distance_to_crystal( - "lehr" - ) - print(f"{pos=}") - print(f"{crystal_dist=}") - print(f"{psd=}") - - arfs = [] - for det in detectors: - # arf actor - arf = sim.add_actor("ARFActor", f"arf_{det.name}") - arf.mother = det.name - arf.batch_size = 1e5 - arf.image_size = [128, 128] - arf.image_spacing = [3 * mm, 3 * mm] - arf.verbose_batch = True - arf.distance_to_crystal = crystal_dist - arf.gpu_mode = "auto" - arf.pth_filename = pth_filename - arf.enable_hit_slice = False - arfs.append(arf) - - return arfs[0], arfs[1] +def add_arf_actor(sim, detector_plane, size, spacing, crystal_dist, name, pth_filename): + # arf actor + arf = sim.add_actor("ARFActor", f"arf_{name}") + arf.mother = detector_plane.name + arf.batch_size = 1e5 + arf.image_size = size + arf.image_spacing = spacing + arf.verbose_batch = True + arf.distance_to_crystal = crystal_dist + arf.gpu_mode = "auto" + arf.enable_hit_slice = False + arf.pth_filename = pth_filename + return arf def set_duration(sim, total_activity, w, duration): From 513b68a85cc85e62f1aa448a8f822b3a6bd9ff57 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 17 Nov 2023 09:21:58 +0100 Subject: [PATCH 09/37] orientation test --- .../test066_spect_gaga_garf_0_orientation.py | 73 +++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100755 opengate/tests/src/test066_spect_gaga_garf_0_orientation.py diff --git a/opengate/tests/src/test066_spect_gaga_garf_0_orientation.py b/opengate/tests/src/test066_spect_gaga_garf_0_orientation.py new file mode 100755 index 000000000..1bb16f7da --- /dev/null +++ b/opengate/tests/src/test066_spect_gaga_garf_0_orientation.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from test066_spect_gaga_garf_helpers import * +from opengate.tests import utility + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, None, "test066") + output_path = paths.output + + # create the simulation + sim = gate.Simulation() + simu_name = "test066_0_orientation" + + # options + ui = sim.user_info + ui.number_of_threads = 1 + ui.visu = True + ui.visu_type = "vrml" + ui.random_seed = "auto" + ui.check_volumes_overlap = True + + # units + mm = gate.g4_units.mm + sec = gate.g4_units.second + Bq = gate.g4_units.Bq + cm3 = gate.g4_units.cm3 + BqmL = Bq / cm3 + + # options + colli_type = "lehr" + radius = 400 * mm + gantry_angle = -20 + + # main elements : spect + phantom + head, crystal = genm670.add_ge_nm67_spect_head( + sim, collimator_type=colli_type, debug=True + ) + + # rotation by default + genm670.set_head_orientation(head, colli_type, radius, gantry_angle) + print(f"Head translation = {head.translation[0]=:.2f}") + + # add arf plane (put in parallel world to avoid overlap) + sim.add_parallel_world("arf_world") + size = [128, 128] + size = [198, 158] # enlarged to see better + spacing = [4 * mm, 4 * mm] + plane_size = [size[0] * spacing[0], size[1] * spacing[1]] + arf_plane = genm670.add_detection_plane_for_arf( + sim, plane_size, colli_type, radius, gantry_angle + ) + arf_plane.mother = "arf_world" + print(f"Plane translation = {arf_plane.translation[0]:.2f}") + + # table + genm670.add_fake_table(sim) + + # add IEC phantom + gate_iec.add_iec_phantom(sim, name="iec") + sim.set_production_cut("iec", "all", 1 * mm) + + # run + sim.run() + + # colors + red = [1, 0.7, 0.7, 0.8] + blue = [0.5, 0.5, 1, 0.8] + gray = [0.5, 0.5, 0.5, 1] + white = [1, 1, 1, 1] + yellow = [1, 1, 0, 1] + green = [0, 1, 0, 1] From 3423280e62061462d7e511fd8c01d062a4b206f8 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 17 Nov 2023 09:24:23 +0100 Subject: [PATCH 10/37] update tol --- opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py b/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py index 8566290bf..e9283372c 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py +++ b/opengate/tests/src/test066_spect_gaga_garf_5_analyse2.py @@ -54,7 +54,7 @@ test_name, axis="x", scaleImageValuesFactor=options.scaling, - sum_tolerance=15, + sum_tolerance=16, tolerance=120, ) and is_ok From 7b809ea983235f2248fcf69b19472faaaa6f23dc Mon Sep 17 00:00:00 2001 From: David Date: Tue, 28 Nov 2023 09:03:56 +0100 Subject: [PATCH 11/37] add arch info in stat actor --- opengate/actors/miscactors.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/opengate/actors/miscactors.py b/opengate/actors/miscactors.py index 78b2d2102..cc2e9501f 100644 --- a/opengate/actors/miscactors.py +++ b/opengate/actors/miscactors.py @@ -4,9 +4,8 @@ import uproot import numpy as np import time - +import platform import opengate_core as g4 - from .base import ActorBase from ..exception import fatal from ..geometry.utility import rot_np_as_g4, vec_np_as_g4, vec_g4_as_np @@ -49,6 +48,7 @@ def __init__(self, user_info=None): self.counts.stop_time = 0 self.counts.init = 0 self.counts.track_types = {} + self.nb_thread = 1 def __del__(self): # print("del SimulationStatisticsActor", self.user_info.name) @@ -75,14 +75,6 @@ def sps(self): return self.counts.step_count / self.counts.duration * sec return 0 - @property - def nb_thread(self): - if self.simulation is not None: - thread = self.simulation.user_info.number_of_threads - else: - thread = "?" - return thread - @property def simu_start_time(self): if not self.simulation is None: @@ -113,16 +105,22 @@ def __str__(self): f"PPS {self.pps:.0f}\n" f"TPS {self.tps:.0f}\n" f"SPS {self.sps:.0f}\n" - f"start {self.counts.start_time}\n" - f"stop {self.counts.stop_time}\n" + f"Start {self.counts.start_time}\n" + f"Stop {self.counts.stop_time}\n" f'Sim start {g4.G4BestUnit(self.simu_start_time, "Time")}\n' f'Sim end {g4.G4BestUnit(self.simu_end_time, "Time")}\n' - f"Threads {self.nb_thread}" + f"Threads {self.nb_thread}\n" + f"Arch {platform.system()}\n" + f"Python {platform.python_version()}" ) if self.user_info.track_types_flag: s += f"\n" f"Track types: {self.counts.track_types}" return s + def StartSimulationAction(self): + g4.GateSimulationStatisticsActor.StartSimulationAction(self) + self.nb_thread = self.simulation.user_info.number_of_threads + def EndSimulationAction(self): g4.GateSimulationStatisticsActor.EndSimulationAction(self) self.counts = Box(self.GetCounts()) @@ -166,6 +164,8 @@ def write(self, filename): s += f"# SPS (Step per sec) = {self.sps:.0f}\n" s += f"# Threads = {self.nb_thread}\n" s += f"# Date = {datetime.now()}\n" + s += f"# Arch = {platform.system()}\n" + s += f"# Python = {platform.python_version()}\n" if self.user_info.track_types_flag: s += f"# Track types:\n" for t in self.counts.track_types: From eaadc77b95cce146a1093c6f0b4ed08679817f45 Mon Sep 17 00:00:00 2001 From: David Date: Tue, 28 Nov 2023 09:05:12 +0100 Subject: [PATCH 12/37] update for parallel world --- opengate/geometry/utility.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/opengate/geometry/utility.py b/opengate/geometry/utility.py index 54c1075e5..6379caeaf 100644 --- a/opengate/geometry/utility.py +++ b/opengate/geometry/utility.py @@ -170,6 +170,7 @@ def get_transform_world_to_local(vol_name): ctr = [0, 0, 0] crot = Rotation.identity().as_matrix() first = True + # FIXME for parallel world while vol_name != __world_name__: pv = g4.G4PhysicalVolumeStore.GetInstance().GetVolume(vol_name, False) tr = vec_g4_as_np(pv.GetObjectTranslation()) @@ -181,7 +182,11 @@ def get_transform_world_to_local(vol_name): else: crot = np.matmul(rot, crot) ctr = rot.dot(ctr) + tr - vol_name = pv.GetMotherLogical().GetName() + + if pv.GetMotherLogical() == None: + vol_name = __world_name__ + else: + vol_name = pv.GetMotherLogical().GetName() return ctr, crot From 9bd5cf5b9ebc761f24bb9e9210e87d357324ccc2 Mon Sep 17 00:00:00 2001 From: David Date: Tue, 28 Nov 2023 09:06:16 +0100 Subject: [PATCH 13/37] add translation param --- opengate/sources/gansources.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opengate/sources/gansources.py b/opengate/sources/gansources.py index c4770081a..56cc7f361 100644 --- a/opengate/sources/gansources.py +++ b/opengate/sources/gansources.py @@ -195,6 +195,7 @@ def __init__( # options self.compute_directions = False self.use_activity_origin = use_activity_origin + self.translation = [0, 0, 0] # variables self.image = None self.cdf_x = self.cdf_y = self.cdf_z = None @@ -242,7 +243,7 @@ def generate_condition(self, n): # x,y,z are in the image coord system # tey are offset according to the coord system (image center or image offset) - p = np.column_stack((x, y, z)) + self.points_offset + p = np.column_stack((x, y, z)) + self.points_offset + self.translation # need direction ? if self.compute_directions is False: From 0b9744188e197bb692ffe4ad5552655acdf0d006 Mon Sep 17 00:00:00 2001 From: David Date: Tue, 28 Nov 2023 09:06:37 +0100 Subject: [PATCH 14/37] correct crop width --- opengate/tests/utility.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/opengate/tests/utility.py b/opengate/tests/utility.py index 3e005e3d7..55186eb24 100644 --- a/opengate/tests/utility.py +++ b/opengate/tests/utility.py @@ -65,6 +65,12 @@ def read_stat_file(filename): stat.counts.track_types = {} if "Date" in line: stat.date = line[len("# Date =") :] + if "Threads" in line: + a = line[len(f"# Threads =") :] + try: + stat.nb_thread = int(a) + except: + stat.nb_thread = "?" return stat @@ -1641,14 +1647,14 @@ def np_plot_slice( def np_plot_slice_h_line(ax, hline, crop_center, crop_width): - x = np.arange(0, 100) + x = np.arange(0, crop_width[0]) c = int(hline - (crop_center[1] - crop_width[1] / 2)) y = [c] * len(x) ax.plot(x, y, color="r") def np_plot_slice_v_line(ax, vline, crop_center, crop_width): - x = np.arange(0, 100) + x = np.arange(0, crop_width[1]) c = int(vline - (crop_center[0] - crop_width[0] / 2)) y = [c] * len(x) ax.plot(y, x, color="r") From 49ff322639c75ff1364048cee61290099b18aae6 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 1 Dec 2023 17:24:02 +0100 Subject: [PATCH 15/37] add arf plane side check --- .../opengate_lib/GateARFActor.cpp | 25 +++++++++++++------ .../opengate_core/opengate_lib/GateARFActor.h | 3 +++ .../opengate_lib/pyGateARFActor.cpp | 3 ++- opengate/actors/arfactors.py | 1 + 4 files changed, 24 insertions(+), 8 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateARFActor.cpp b/core/opengate_core/opengate_lib/GateARFActor.cpp index 5bceae2f5..ad9ec34ee 100644 --- a/core/opengate_core/opengate_lib/GateARFActor.cpp +++ b/core/opengate_core/opengate_lib/GateARFActor.cpp @@ -36,6 +36,7 @@ void GateARFActor::EndOfRunAction(const G4Run * /*run*/) { l.fPositionY.clear(); l.fDirectionX.clear(); l.fDirectionY.clear(); + l.fDirectionZ.clear(); l.fCurrentNumberOfHits = 0; } } @@ -43,10 +44,21 @@ void GateARFActor::EndOfRunAction(const G4Run * /*run*/) { void GateARFActor::SteppingAction(G4Step *step) { auto &l = fThreadLocalData.Get(); + // get direction and transform to local + auto *pre = step->GetPreStepPoint(); + auto dir = pre->GetMomentumDirection(); + dir = pre->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(dir); + + // which side of the plane ? + if (dir[2] < 0) + return; + l.fCurrentNumberOfHits++; + l.fDirectionX.push_back(dir[0]); + l.fDirectionY.push_back(dir[1]); + // l.fDirectionZ.push_back(dir[2]); // not used // get energy - auto *pre = step->GetPreStepPoint(); l.fEnergy.push_back(pre->GetKineticEnergy()); // get position and transform to local @@ -56,12 +68,6 @@ void GateARFActor::SteppingAction(G4Step *step) { l.fPositionX.push_back(pos[0]); l.fPositionY.push_back(pos[1]); - // get direction and transform to local - auto dir = pre->GetMomentumDirection(); - dir = pre->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(dir); - l.fDirectionX.push_back(dir[0]); - l.fDirectionY.push_back(dir[1]); - // trigger the "apply" (ARF) if the number of hits in the batch is reached if (l.fCurrentNumberOfHits >= fBatchSize) { fApply(this); @@ -70,6 +76,7 @@ void GateARFActor::SteppingAction(G4Step *step) { l.fPositionY.clear(); l.fDirectionX.clear(); l.fDirectionY.clear(); + l.fDirectionZ.clear(); l.fCurrentNumberOfHits = 0; } } @@ -101,3 +108,7 @@ std::vector GateARFActor::GetDirectionX() const { std::vector GateARFActor::GetDirectionY() const { return fThreadLocalData.Get().fDirectionY; } + +std::vector GateARFActor::GetDirectionZ() const { + return fThreadLocalData.Get().fDirectionZ; +} diff --git a/core/opengate_core/opengate_lib/GateARFActor.h b/core/opengate_core/opengate_lib/GateARFActor.h index 306d74afc..7f0902858 100644 --- a/core/opengate_core/opengate_lib/GateARFActor.h +++ b/core/opengate_core/opengate_lib/GateARFActor.h @@ -43,6 +43,8 @@ class GateARFActor : public GateVActor { std::vector GetDirectionY() const; + std::vector GetDirectionZ() const; + // Main function called every step in attached volume void SteppingAction(G4Step *) override; @@ -60,6 +62,7 @@ class GateARFActor : public GateVActor { std::vector fPositionY; std::vector fDirectionX; std::vector fDirectionY; + std::vector fDirectionZ; // number of particle hitting the detector int fCurrentNumberOfHits; // Current run id (to detect if run has changed) diff --git a/core/opengate_core/opengate_lib/pyGateARFActor.cpp b/core/opengate_core/opengate_lib/pyGateARFActor.cpp index 6683d6b38..80949c1cc 100644 --- a/core/opengate_core/opengate_lib/pyGateARFActor.cpp +++ b/core/opengate_core/opengate_lib/pyGateARFActor.cpp @@ -24,5 +24,6 @@ void init_GateARFActor(py::module &m) { .def("GetPositionX", &GateARFActor::GetPositionX) .def("GetPositionY", &GateARFActor::GetPositionY) .def("GetDirectionX", &GateARFActor::GetDirectionX) - .def("GetDirectionY", &GateARFActor::GetDirectionY); + .def("GetDirectionY", &GateARFActor::GetDirectionY) + .def("GetDirectionZ", &GateARFActor::GetDirectionZ); } diff --git a/opengate/actors/arfactors.py b/opengate/actors/arfactors.py index edd2f1e85..140bf0f48 100644 --- a/opengate/actors/arfactors.py +++ b/opengate/actors/arfactors.py @@ -239,6 +239,7 @@ def apply_with_lock(self, actor): print( f"Apply ARF to {energy.shape[0]} hits (device = {self.model_data['current_gpu_mode']})" ) + ax = x[:, 2:5] # two angles and energy # FIXME index ? w = self.garf.nn_predict(self.model, self.nn["model_data"], ax) From 9c0f8ce86064a309673282d76616dca2c4328e45 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 1 Dec 2023 17:25:27 +0100 Subject: [PATCH 16/37] add src translation --- core/opengate_core/opengate_lib/GateGANSource.cpp | 1 + opengate/sources/gansources.py | 6 ++++++ 2 files changed, 7 insertions(+) diff --git a/core/opengate_core/opengate_lib/GateGANSource.cpp b/core/opengate_core/opengate_lib/GateGANSource.cpp index 8e4e193b2..949d8aa9a 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANSource.cpp @@ -224,6 +224,7 @@ G4ThreeVector GateGANSource::GeneratePrimariesPosition() { position = G4ThreeVector(fPositionX[fCurrentIndex], fPositionY[fCurrentIndex], fPositionZ[fCurrentIndex]); + position = fLocalRotation * position + fLocalTranslation; // FIXME // move position according to mother volume auto &l = fThreadLocalData.Get(); position = l.fGlobalRotation * position + l.fGlobalTranslation; diff --git a/opengate/sources/gansources.py b/opengate/sources/gansources.py index 56cc7f361..1f82a1733 100644 --- a/opengate/sources/gansources.py +++ b/opengate/sources/gansources.py @@ -12,6 +12,7 @@ from ..image import get_info_from_image from ..image import compute_image_3D_CDF from .generic import generate_isotropic_directions +from scipy.spatial.transform import Rotation def import_gaga_phsp(): @@ -196,6 +197,7 @@ def __init__( self.compute_directions = False self.use_activity_origin = use_activity_origin self.translation = [0, 0, 0] + self.rotation = Rotation.identity().as_matrix() # variables self.image = None self.cdf_x = self.cdf_y = self.cdf_z = None @@ -245,10 +247,14 @@ def generate_condition(self, n): # tey are offset according to the coord system (image center or image offset) p = np.column_stack((x, y, z)) + self.points_offset + self.translation + # rotation + p = np.dot(p, self.rotation.T) + # need direction ? if self.compute_directions is False: return p v = generate_isotropic_directions(n, rs=self.rs) + v = np.dot(v, self.rotation.T) return np.column_stack((p, v)) From 013015b5f6ac9fd012d4d3758262902210e52a80 Mon Sep 17 00:00:00 2001 From: David Date: Mon, 4 Dec 2023 11:12:40 +0100 Subject: [PATCH 17/37] correct apply g4 command --- core/external/pybind11 | 2 +- opengate/engines.py | 11 ++- opengate/geometry/utility.py | 1 - opengate/managers.py | 2 +- opengate/tests/src/test004_simple.py | 5 +- opengate/tests/src/test004_simple_mt.py | 4 +- .../tests/src/test004_simple_visu_gdml.py | 2 +- opengate/tests/src/test004_simple_visu_qt.py | 2 +- .../src/test004_simple_visu_vrml_file.py | 2 +- opengate/tests/src/test005_proton.py | 6 +- opengate/tests/src/test010_generic_source.py | 6 +- opengate/tests/src/test012_mt_dose_actor.py | 6 +- opengate/tests/src/test013_phys_lists_1.py | 2 +- .../tests/src/test013_phys_lists_2_wip.py | 2 +- .../tests/src/test013_phys_lists_3_wip.py | 2 +- opengate/tests/src/test013_phys_lists_4.py | 2 +- .../tests/src/test013_phys_lists_5_wip.py | 2 +- .../tests/src/test015_iec_phantom_4_wip.py | 93 +++++++++++++++++++ opengate/tests/src/test026_simple_signal.py | 2 +- .../test066_spect_gaga_garf_0_orientation.py | 2 +- .../test066_spect_gaga_garf_1_reference.py | 2 +- .../src/test066_spect_gaga_garf_helpers.py | 6 +- 22 files changed, 130 insertions(+), 34 deletions(-) create mode 100755 opengate/tests/src/test015_iec_phantom_4_wip.py diff --git a/core/external/pybind11 b/core/external/pybind11 index a34596bfe..a67d78657 160000 --- a/core/external/pybind11 +++ b/core/external/pybind11 @@ -1 +1 @@ -Subproject commit a34596bfe1947b4a6b0bcc4218e1f72d0c2e9b4c +Subproject commit a67d786571cc6d2c2ecec33985e00e6d22dc68e5 diff --git a/opengate/engines.py b/opengate/engines.py index 472f2853a..2460a31e7 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -1472,7 +1472,7 @@ def initialize_random_engine(self): # set the seed g4.G4Random.setTheSeed(self.current_random_seed, 0) - def initialize_g4_verbose(self): + def initialize_g4_verbose(self, verbose_tracking=True): if self.simulation.user_info.g4_verbose: # Geant4 output with color ui = UIsessionVerbose() @@ -1480,9 +1480,12 @@ def initialize_g4_verbose(self): else: # no Geant4 output ui = UIsessionSilent() - self.simulation.add_g4_command_after_init( - f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" - ) + + # do not set verbose tracking when "overlap" + if self.simulation.g4_verbose_level_tracking >= 0: + self.simulation.add_g4_command_after_init( + f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" + ) # it is also possible to set ui=None for 'default' output # we must keep a ref to ui_session self.ui_session = ui diff --git a/opengate/geometry/utility.py b/opengate/geometry/utility.py index 863b7563a..a1561586e 100644 --- a/opengate/geometry/utility.py +++ b/opengate/geometry/utility.py @@ -188,7 +188,6 @@ def get_transform_world_to_local(volume): ctr = volume.translation_list[i] crot = volume.rotation_list[i] for vol in volume.ancestor_volumes[::-1]: - print(vol.name) crot = np.matmul(vol.rotation_list[0], crot) ctr = vol.rotation_list[0].dot(ctr) + vol.translation_list[0] cumulative_translation.append(ctr) diff --git a/opengate/managers.py b/opengate/managers.py index 4231786d9..3a15c5afc 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -909,7 +909,7 @@ class Simulation(GateObject): ), "g4_verbose": (False, {"doc": "Switch on/off Geant4's verbose output."}), "g4_verbose_level_tracking": ( - 0, + -1, { "doc": "Activate verbose tracking in Geant4 via G4 command '/tracking/verbose g4_verbose_level_tracking'." }, diff --git a/opengate/tests/src/test004_simple.py b/opengate/tests/src/test004_simple.py index 6d60afc4c..2b4be98a0 100755 --- a/opengate/tests/src/test004_simple.py +++ b/opengate/tests/src/test004_simple.py @@ -109,9 +109,10 @@ """ Start the simulation ! You can relax and drink coffee. - (The commented line indicates how to indicate to Geant4 to verbose during the simulation). + (The commented line indicates how to indicate to Geant4 to verbose during the simulation, + if the flag sim.g4_verbose is True). """ - # sim.apply_g4_command("/run/verbose 1") + # sim.add_g4_command_after_init("/run/verbose 1") sim.user_fct_after_init = check_production_cuts sim.run() diff --git a/opengate/tests/src/test004_simple_mt.py b/opengate/tests/src/test004_simple_mt.py index 9aac6a511..9e38b8945 100755 --- a/opengate/tests/src/test004_simple_mt.py +++ b/opengate/tests/src/test004_simple_mt.py @@ -66,8 +66,8 @@ s.track_types_flag = True # start simulation - # sim.apply_g4_command("/run/verbose 0") - # sim.apply_g4_command("/run/eventModulo 5000 1") + # sim.add_g4_command_after_init("/run/verbose 0") + # sim.add_g4_command_after_init("/run/eventModulo 5000 1") sim.run() # get results diff --git a/opengate/tests/src/test004_simple_visu_gdml.py b/opengate/tests/src/test004_simple_visu_gdml.py index 4232d99f1..9effdbfd2 100755 --- a/opengate/tests/src/test004_simple_visu_gdml.py +++ b/opengate/tests/src/test004_simple_visu_gdml.py @@ -55,7 +55,7 @@ sim.add_actor("SimulationStatisticsActor", "Stats") # start simulation - # sim.apply_g4_command("/run/verbose 1") + # sim.add_g4_command_after_init("/run/verbose 1") sim.run() stats = sim.output.get_actor("Stats") diff --git a/opengate/tests/src/test004_simple_visu_qt.py b/opengate/tests/src/test004_simple_visu_qt.py index 37c8c70a3..f4ecd6f42 100755 --- a/opengate/tests/src/test004_simple_visu_qt.py +++ b/opengate/tests/src/test004_simple_visu_qt.py @@ -53,7 +53,7 @@ sim.add_actor("SimulationStatisticsActor", "Stats") # start simulation - # sim.apply_g4_command("/run/verbose 1") + # sim.add_g4_command_after_init("/run/verbose 1") sim.run() stats = sim.output.get_actor("Stats") diff --git a/opengate/tests/src/test004_simple_visu_vrml_file.py b/opengate/tests/src/test004_simple_visu_vrml_file.py index 35e2ab9ba..be890facc 100755 --- a/opengate/tests/src/test004_simple_visu_vrml_file.py +++ b/opengate/tests/src/test004_simple_visu_vrml_file.py @@ -56,7 +56,7 @@ sim.add_actor("SimulationStatisticsActor", "Stats") # start simulation - # sim.apply_g4_command("/run/verbose 1") + # sim.add_g4_command_after_init("/run/verbose 1") sim.run(True) # visu diff --git a/opengate/tests/src/test005_proton.py b/opengate/tests/src/test005_proton.py index 73006ce6a..e9a161a41 100755 --- a/opengate/tests/src/test005_proton.py +++ b/opengate/tests/src/test005_proton.py @@ -41,9 +41,9 @@ # verbose (WARNING : sim.g4_verbose must be True !) sim.add_g4_command_after_init("/tracking/verbose 0") - # sim.apply_g4_command("/run/verbose 2") - # sim.apply_g4_command("/event/verbose 2") - # sim.apply_g4_command("/tracking/verbose 1") + # sim.add_g4_command_after_init("/run/verbose 2") + # sim.add_g4_command_after_init("/event/verbose 2") + # sim.add_g4_command_after_init("/tracking/verbose 1") print(sim.source_manager.dump_sources()) diff --git a/opengate/tests/src/test010_generic_source.py b/opengate/tests/src/test010_generic_source.py index cc0b897c9..42f00ecc9 100755 --- a/opengate/tests/src/test010_generic_source.py +++ b/opengate/tests/src/test010_generic_source.py @@ -101,9 +101,9 @@ # verbose sim.add_g4_command_after_init("/tracking/verbose 0") - # sim.apply_g4_command("/run/verbose 2") - # sim.apply_g4_command("/event/verbose 2") - # sim.apply_g4_command("/tracking/verbose 1") + # sim.add_g4_command_after_init("/run/verbose 2") + # sim.add_g4_command_after_init("/event/verbose 2") + # sim.add_g4_command_after_init("/tracking/verbose 1") # start simulation sim.run() diff --git a/opengate/tests/src/test012_mt_dose_actor.py b/opengate/tests/src/test012_mt_dose_actor.py index 3eaa760b5..ab3e376bc 100755 --- a/opengate/tests/src/test012_mt_dose_actor.py +++ b/opengate/tests/src/test012_mt_dose_actor.py @@ -80,10 +80,10 @@ print(sim.volume_manager.dump_volumes()) # verbose - # sim.apply_g4_command('/tracking/verbose 0') + # sim.add_g4_command_after_init('/tracking/verbose 0') sim.add_g4_command_after_init("/run/verbose 2") - # sim.apply_g4_command("/event/verbose 2") - # sim.apply_g4_command("/tracking/verbose 1") + # sim.add_g4_command_after_init("/event/verbose 2") + # sim.add_g4_command_after_init("/tracking/verbose 1") # start simulation sim.run() diff --git a/opengate/tests/src/test013_phys_lists_1.py b/opengate/tests/src/test013_phys_lists_1.py index ff837f52b..01182ea99 100755 --- a/opengate/tests/src/test013_phys_lists_1.py +++ b/opengate/tests/src/test013_phys_lists_1.py @@ -36,7 +36,7 @@ print(sim.physics_manager.dump_production_cuts()) # start simulation - # sim.apply_g4_command("/tracking/verbose 1") + # sim.add_g4_command_after_init("/tracking/verbose 1") sim.run() # Gate mac/main_1.mac diff --git a/opengate/tests/src/test013_phys_lists_2_wip.py b/opengate/tests/src/test013_phys_lists_2_wip.py index 68f216b2d..5978a250d 100755 --- a/opengate/tests/src/test013_phys_lists_2_wip.py +++ b/opengate/tests/src/test013_phys_lists_2_wip.py @@ -32,7 +32,7 @@ # start simulation # sim.set_g4_verbose(True) -# sim.apply_g4_command("/tracking/verbose 1") +# sim.add_g4_command_after_init("/tracking/verbose 1") sim.run() stats = sim.output.get_actor("Stats") diff --git a/opengate/tests/src/test013_phys_lists_3_wip.py b/opengate/tests/src/test013_phys_lists_3_wip.py index 51e5adae9..29c0052b0 100755 --- a/opengate/tests/src/test013_phys_lists_3_wip.py +++ b/opengate/tests/src/test013_phys_lists_3_wip.py @@ -44,7 +44,7 @@ # start simulation sim.g4_verbose = False - # sim.apply_g4_command("/tracking/verbose 1") + # sim.add_g4_command_after_init("/tracking/verbose 1") sim.user_fct_after_init = check_production_cuts sim.run() diff --git a/opengate/tests/src/test013_phys_lists_4.py b/opengate/tests/src/test013_phys_lists_4.py index 7931aabd7..8fa8fca47 100755 --- a/opengate/tests/src/test013_phys_lists_4.py +++ b/opengate/tests/src/test013_phys_lists_4.py @@ -42,7 +42,7 @@ # start simulation sim.g4_verbose = True - # sim.apply_g4_command("/tracking/verbose 1") + # sim.add_g4_command_after_init("/tracking/verbose 1") sim.run() stats = sim.output.get_actor("Stats") diff --git a/opengate/tests/src/test013_phys_lists_5_wip.py b/opengate/tests/src/test013_phys_lists_5_wip.py index 2409e2875..11026c363 100755 --- a/opengate/tests/src/test013_phys_lists_5_wip.py +++ b/opengate/tests/src/test013_phys_lists_5_wip.py @@ -51,7 +51,7 @@ print(sim.physics_manager.dump_production_cuts()) # start simulation -# sim.apply_g4_command("/tracking/verbose 1") +# sim.add_g4_command_after_init("/tracking/verbose 1") sim.g4_verbose = False sim.user_fct_after_init = check_production_cuts sim.run() diff --git a/opengate/tests/src/test015_iec_phantom_4_wip.py b/opengate/tests/src/test015_iec_phantom_4_wip.py new file mode 100755 index 000000000..15419fc84 --- /dev/null +++ b/opengate/tests/src/test015_iec_phantom_4_wip.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +from scipy.spatial.transform import Rotation +import opengate.contrib.phantoms.nemaiec as gate_iec + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "", "test015") + + # create the simulation + sim = gate.Simulation() + + # units + MeV = gate.g4_units.MeV + m = gate.g4_units.m + mm = gate.g4_units.mm + cm = gate.g4_units.cm + cm3 = gate.g4_units.cm3 + Bq = gate.g4_units.Bq + BqmL = Bq / cm3 + + # main options + ui = sim.user_info + # ui.visu = True + # ui.visu_type = "qt" + ui.check_volumes_overlap = True + ui.random_seed = 123654987 + + # physics + p = sim.get_physics_user_info() + p.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.set_production_cut("world", "all", 10 * mm) + + # world size + world = sim.world + world.size = [0.5 * m, 0.5 * m, 0.5 * m] + + # add an iec phantom + # rotation 180 around X to be like in the iec 61217 coordinate system + iec_phantom = gate_iec.add_iec_phantom(sim) + + # box around + b = sim.add_volume("Box", "box") + b.size = [33 * cm, 33 * cm, 25 * cm] + iec_phantom.mother = b.name + + # add sources for all spheres + a = 100 * BqmL + activity_Bq_mL = [10 * a, 2 * a, 3 * a, 4 * a, 5 * a, 6 * a] + sources = gate_iec.add_spheres_sources( + sim, iec_phantom.name, "sources", "all", activity_Bq_mL, verbose=True + ) + for source in sources: + source.particle = "alpha" + source.energy.type = "mono" + source.energy.mono = 100 * MeV + + # add stat actor + stats = sim.add_actor("SimulationStatisticsActor", "stats") + stats.track_types_flag = True + stats.output = paths.output / "test015_iec_2_stats.txt" + + # add dose actor + dose = sim.add_actor("DoseActor", "dose") + dose.output = paths.output / "test015_iec_2.mhd" + dose.mother = iec_phantom.name + dose.size = [100, 100, 100] + dose.spacing = [2 * mm, 2 * mm, 2 * mm] + + # start + sim.run() + + # compare stats + stats = sim.output.get_actor("stats") + stats_ref = utility.read_stat_file(paths.output_ref / "test015_iec_2_stats.txt") + is_ok = utility.assert_stats(stats, stats_ref, tolerance=0.03) + + # compare images + f = paths.output / "test015_iec_2.mhd" + im_ok = utility.assert_images( + paths.output_ref / "test015_iec_2.mhd", + dose.output, + stats, + axis="x", + tolerance=40, + ignore_value=0, + sum_tolerance=2, + ) + + is_ok = is_ok and im_ok + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test026_simple_signal.py b/opengate/tests/src/test026_simple_signal.py index 7f4fbf5e0..65ce60a20 100755 --- a/opengate/tests/src/test026_simple_signal.py +++ b/opengate/tests/src/test026_simple_signal.py @@ -48,7 +48,7 @@ stats.track_types_flag = True # start simulation - # sim.apply_g4_command("/run/verbose 1") + # sim.add_g4_command_after_init("/run/verbose 1") sim.run() # get result diff --git a/opengate/tests/src/test066_spect_gaga_garf_0_orientation.py b/opengate/tests/src/test066_spect_gaga_garf_0_orientation.py index 1bb16f7da..80ec47037 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_0_orientation.py +++ b/opengate/tests/src/test066_spect_gaga_garf_0_orientation.py @@ -59,7 +59,7 @@ # add IEC phantom gate_iec.add_iec_phantom(sim, name="iec") - sim.set_production_cut("iec", "all", 1 * mm) + sim.physics_manager.set_production_cut("iec", "all", 1 * mm) # run sim.run() diff --git a/opengate/tests/src/test066_spect_gaga_garf_1_reference.py b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py index 836e70e7d..4a3cccb20 100755 --- a/opengate/tests/src/test066_spect_gaga_garf_1_reference.py +++ b/opengate/tests/src/test066_spect_gaga_garf_1_reference.py @@ -52,7 +52,7 @@ # add IEC phantom iec = gate_iec.add_iec_phantom(sim, name="iec") - sim.set_production_cut("iec", "all", 1 * mm) + sim.physics_manager.set_production_cut("iec", "all", 1 * mm) iec.rotation = Rotation.from_euler("x", 90, degrees=True).as_matrix() # sources IEC diff --git a/opengate/tests/src/test066_spect_gaga_garf_helpers.py b/opengate/tests/src/test066_spect_gaga_garf_helpers.py index dd2e7f557..22cc1ef61 100644 --- a/opengate/tests/src/test066_spect_gaga_garf_helpers.py +++ b/opengate/tests/src/test066_spect_gaga_garf_helpers.py @@ -20,7 +20,7 @@ def set_phys(sim): m = gate.g4_units.m p = sim.get_physics_user_info() p.physics_list_name = "G4EmStandardPhysics_option3" - sim.set_production_cut("world", "all", 1e3 * m) + sim.physics_manager.set_production_cut("world", "all", 1e3 * m) def create_simu_with_genm670(sim, collimator_type="lehr", debug=False): @@ -59,8 +59,8 @@ def create_simu_with_genm670(sim, collimator_type="lehr", debug=False): # phys set_phys(sim) - sim.set_production_cut("spect1", "all", 1 * mm) - sim.set_production_cut("spect2", "all", 1 * mm) + sim.physics_manager.set_production_cut("spect1", "all", 1 * mm) + sim.physics_manager.set_production_cut("spect2", "all", 1 * mm) return proj1, proj2 From 57f4e323ebae029c0e5af15f797621dfede4b417 Mon Sep 17 00:00:00 2001 From: David Date: Mon, 4 Dec 2023 11:48:55 +0100 Subject: [PATCH 18/37] remove unused param --- opengate/engines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/engines.py b/opengate/engines.py index 2460a31e7..12c14bbbc 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -1472,7 +1472,7 @@ def initialize_random_engine(self): # set the seed g4.G4Random.setTheSeed(self.current_random_seed, 0) - def initialize_g4_verbose(self, verbose_tracking=True): + def initialize_g4_verbose(self): if self.simulation.user_info.g4_verbose: # Geant4 output with color ui = UIsessionVerbose() From b18c46cb4e737e4e0f311fc29f017a83de5dc01a Mon Sep 17 00:00:00 2001 From: David Date: Mon, 4 Dec 2023 11:49:30 +0100 Subject: [PATCH 19/37] update fct name --- opengate/tests/src/test034_gan_phsp_linac.py | 2 +- .../src/test038_gan_phsp_spect_gan_helpers.py | 2 +- opengate/tests/src/test040_gan_phsp_pet_gan.py | 2 +- opengate/tests/src/test043_garf.py | 18 ++++++++++-------- opengate/tests/src/test043_garf_helpers.py | 2 +- opengate/tests/src/test043_garf_mt.py | 2 +- .../src/test045_gan_phsp_pet_gan_helpers.py | 2 +- .../tests/src/test047_gan_vox_source_cond.py | 2 +- opengate/tests/utility.py | 5 +++-- 9 files changed, 20 insertions(+), 17 deletions(-) diff --git a/opengate/tests/src/test034_gan_phsp_linac.py b/opengate/tests/src/test034_gan_phsp_linac.py index 9c8231661..f93bbbb39 100755 --- a/opengate/tests/src/test034_gan_phsp_linac.py +++ b/opengate/tests/src/test034_gan_phsp_linac.py @@ -73,7 +73,7 @@ # it is possible to define another generator # gsource.generator = generator gsource.gpu_mode = ( - utility.get_gpu_mode() + utility.get_gpu_mode_for_tests() ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors # add stat actor diff --git a/opengate/tests/src/test038_gan_phsp_spect_gan_helpers.py b/opengate/tests/src/test038_gan_phsp_spect_gan_helpers.py index 8baaa2379..734649adb 100644 --- a/opengate/tests/src/test038_gan_phsp_spect_gan_helpers.py +++ b/opengate/tests/src/test038_gan_phsp_spect_gan_helpers.py @@ -146,7 +146,7 @@ def create_simulation(sim, paths, colli="lehr"): gsource.batch_size = 5e4 gsource.verbose_generator = True gsource.gpu_mode = ( - utility.get_gpu_mode() + utility.get_gpu_mode_for_tests() ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors # GANSourceConditionalGenerator manages the conditional GAN diff --git a/opengate/tests/src/test040_gan_phsp_pet_gan.py b/opengate/tests/src/test040_gan_phsp_pet_gan.py index 03d4ae15f..20d062252 100755 --- a/opengate/tests/src/test040_gan_phsp_pet_gan.py +++ b/opengate/tests/src/test040_gan_phsp_pet_gan.py @@ -162,7 +162,7 @@ def gen_cond(n): gsource, 210 * mm, gen_cond ) gsource.gpu_mode = ( - utility.get_gpu_mode() + utility.get_gpu_mode_for_tests() ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors # add stat actor diff --git a/opengate/tests/src/test043_garf.py b/opengate/tests/src/test043_garf.py index 92c84ae34..2fad11a53 100755 --- a/opengate/tests/src/test043_garf.py +++ b/opengate/tests/src/test043_garf.py @@ -15,7 +15,8 @@ sim.g4_verbose = False sim.g4_verbose_level = 1 sim.number_of_threads = 1 - sim.visu = False + # sim.visu = True + sim.visu_type = "vrml" sim.random_seed = 321654987 # units @@ -28,6 +29,7 @@ # activity activity = 1e6 * Bq / sim.number_of_threads + # activity = 1e2 * Bq / sim.number_of_threads # add a material database sim.volume_manager.add_material_database( @@ -65,18 +67,18 @@ arf.image_spacing = [4.41806 * mm, 4.41806 * mm] arf.verbose_batch = True arf.distance_to_crystal = crystal_dist # 74.625 * mm - arf.distance_to_crystal = 74.625 * mm + arf.distance_to_crystal = 74.625 * mm ## FIXME arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" arf.enable_hit_slice = True - arf.gpu_mode = ( - utility.get_gpu_mode() - ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors - + # the gpu_mode should be "auto" but when running on github actions, + # with set "cpu" to avoid mps errors + arf.gpu_mode = utility.get_gpu_mode_for_tests() # add stat actor s = sim.add_actor("SimulationStatisticsActor", "stats") s.track_types_flag = True - # start simulation to check if ok with start_new_process + # FIXME + """# start simulation to check if ok with start_new_process s1a = s1.activity s2a = s2.activity s3a = s3.activity @@ -88,7 +90,7 @@ # restart simulation s1.activity = s1a s2.activity = s2a - s3.activity = s3a + s3.activity = s3a""" sim.run(False) # print results at the end diff --git a/opengate/tests/src/test043_garf_helpers.py b/opengate/tests/src/test043_garf_helpers.py index 4bcff4ad1..6bc16ad14 100644 --- a/opengate/tests/src/test043_garf_helpers.py +++ b/opengate/tests/src/test043_garf_helpers.py @@ -140,7 +140,7 @@ def create_sim_test_region(sim): arf.pth_filename = paths.gate_data / "pth" / "arf_Tc99m_v3.pth" arf.enable_hit_slice = True arf.gpu_mode = ( - utility.get_gpu_mode() + utility.get_gpu_mode_for_tests() ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors # add stat actor diff --git a/opengate/tests/src/test043_garf_mt.py b/opengate/tests/src/test043_garf_mt.py index a9deafb4e..fed7555d6 100755 --- a/opengate/tests/src/test043_garf_mt.py +++ b/opengate/tests/src/test043_garf_mt.py @@ -68,7 +68,7 @@ arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" arf.enable_hit_slice = True arf.gpu_mode = ( - utility.get_gpu_mode() + utility.get_gpu_mode_for_tests() ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors # add stat actor diff --git a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py index cd75176a5..01acc527c 100644 --- a/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py +++ b/opengate/tests/src/test045_gan_phsp_pet_gan_helpers.py @@ -245,7 +245,7 @@ def gen_cond(n): gsource, 210 * mm, gen_cond ) gsource.generator = gen - gsource.gpu_mode = utility.get_gpu_mode() + gsource.gpu_mode = utility.get_gpu_mode_for_tests() def add_gaga_source_vox_condition(sim, p): diff --git a/opengate/tests/src/test047_gan_vox_source_cond.py b/opengate/tests/src/test047_gan_vox_source_cond.py index c3387bfa9..c26c978ab 100755 --- a/opengate/tests/src/test047_gan_vox_source_cond.py +++ b/opengate/tests/src/test047_gan_vox_source_cond.py @@ -97,7 +97,7 @@ source.batch_size = 1e5 source.verbose_generator = True source.gpu_mode = ( - utility.get_gpu_mode() + utility.get_gpu_mode_for_tests() ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors # cuts (not need precision here) diff --git a/opengate/tests/utility.py b/opengate/tests/utility.py index 8d027dd82..f6a3f1698 100644 --- a/opengate/tests/utility.py +++ b/opengate/tests/utility.py @@ -1556,10 +1556,11 @@ def compare_trees4(p1, p2, param): return is_ok -def get_gpu_mode(): +def get_gpu_mode_for_tests(): """ return "auto" except if the test runs with macos and github actions - On macos and github actions, mps is detected but not usable and lead to errors. So choose "cpu" in such a case + On macos and github actions, mps is detected but not usable and lead to errors. + So we choose "cpu" in such a case """ if "GITHUB_WORKSPACE" in os.environ and sys.platform == "darwin": print("Detection of Github actions and MacOS -> Use CPU") From 995e3cfdf001965df311924d778757e2720ba158 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 17:55:08 +0100 Subject: [PATCH 20/37] first version --- core/external/pybind11 | 2 +- core/opengate_core/opengate_core.cpp | 3 + .../opengate_lib/GateFluenceActor.cpp | 91 ++++++++++++++ .../opengate_lib/GateFluenceActor.h | 58 +++++++++ .../opengate_lib/GateGANSource.cpp | 1 + .../opengate_lib/pyGateFluenceActor.cpp | 22 ++++ opengate/actors/actorbuilders.py | 3 +- opengate/actors/doseactors.py | 111 +++++++++++++++++- opengate/tests/data | 2 +- 9 files changed, 289 insertions(+), 4 deletions(-) create mode 100644 core/opengate_core/opengate_lib/GateFluenceActor.cpp create mode 100644 core/opengate_core/opengate_lib/GateFluenceActor.h create mode 100644 core/opengate_core/opengate_lib/pyGateFluenceActor.cpp diff --git a/core/external/pybind11 b/core/external/pybind11 index a34596bfe..a67d78657 160000 --- a/core/external/pybind11 +++ b/core/external/pybind11 @@ -1 +1 @@ -Subproject commit a34596bfe1947b4a6b0bcc4218e1f72d0c2e9b4c +Subproject commit a67d786571cc6d2c2ecec33985e00e6d22dc68e5 diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index 23d4a137c..86e314219 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -249,6 +249,8 @@ void init_GateKineticEnergyFilter(py::module &); void init_GateDoseActor(py::module &m); +void init_GateFluenceActor(py::module &m); + void init_GateLETActor(py::module &m); void init_GateARFActor(py::module &m); @@ -475,6 +477,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateEventAction(m); init_GateTrackingAction(m); init_GateDoseActor(m); + init_GateFluenceActor(m); init_GateLETActor(m); init_GateSimulationStatisticsActor(m); init_GatePhaseSpaceActor(m); diff --git a/core/opengate_core/opengate_lib/GateFluenceActor.cpp b/core/opengate_core/opengate_lib/GateFluenceActor.cpp new file mode 100644 index 000000000..bc1938646 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateFluenceActor.cpp @@ -0,0 +1,91 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ------------------------------------ -------------- */ + +#include "G4Navigator.hh" +#include "G4RunManager.hh" +#include "G4Threading.hh" + +#include "GateFluenceActor.h" +#include "GateHelpersDict.h" +#include "GateHelpersImage.h" + +#include +#include +#include + +// Mutex that will be used by thread to write the output image +G4Mutex SetPixelFluenceMutex = G4MUTEX_INITIALIZER; + +GateFluenceActor::GateFluenceActor(py::dict &user_info) + : GateVActor(user_info, true) { + // Create the image pointer + // (the size and allocation will be performed on the py side) + cpp_fluence_image = Image3DType::New(); + // Action for this actor: during stepping + fActions.insert("SteppingAction"); + fActions.insert("BeginOfRunAction"); + fActions.insert("BeginOfEventAction"); + fActions.insert("EndSimulationAction"); + // translation + fInitialTranslation = DictGetG4ThreeVector(user_info, "translation"); +} + +void GateFluenceActor::ActorInitialize() {} + +void GateFluenceActor::BeginOfRunAction(const G4Run *) { + Image3DType::RegionType region = + cpp_fluence_image->GetLargestPossibleRegion(); + size_fluence = region.GetSize(); + + // Important ! The volume may have moved, so we re-attach each run + AttachImageToVolume(cpp_fluence_image, fPhysicalVolumeName, + fInitialTranslation); + // compute volume of a dose voxel + auto sp = cpp_fluence_image->GetSpacing(); + fVoxelVolume = sp[0] * sp[1] * sp[2]; +} + +void GateFluenceActor::BeginOfEventAction(const G4Event *event) {} + +void GateFluenceActor::SteppingAction(G4Step *step) { + // same method to consider only entering tracks + if (step->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) { + // the pre-position is at the edge + auto preGlobal = step->GetPreStepPoint()->GetPosition(); + auto dir = step->GetPreStepPoint()->GetMomentumDirection(); + auto touchable = step->GetPreStepPoint()->GetTouchable(); + + // consider position in the local volume, slightly shifted by 0.1 nm because + // otherwise, it can be considered as outside the volume by isInside. + auto position = preGlobal + 0.1 * CLHEP::nm * dir; + auto localPosition = + touchable->GetHistory()->GetTransform(0).TransformPoint(position); + + // convert G4ThreeVector to itk PointType + Image3DType::PointType point; + point[0] = localPosition[0]; + point[1] = localPosition[1]; + point[2] = localPosition[2]; + + // get weight + auto w = step->GetTrack()->GetWeight(); + + // get pixel index + Image3DType::IndexType index; + bool isInside = + cpp_fluence_image->TransformPhysicalPointToIndex(point, index); + + // set value + if (isInside) { + G4AutoLock FluenceMutex(&SetPixelFluenceMutex); + // add hit + ImageAddValue(cpp_fluence_image, index, w); + } // else : outside the image + } +} + +void GateFluenceActor::EndSimulationAction() {} diff --git a/core/opengate_core/opengate_lib/GateFluenceActor.h b/core/opengate_core/opengate_lib/GateFluenceActor.h new file mode 100644 index 000000000..d296b6d89 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateFluenceActor.h @@ -0,0 +1,58 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateFluenceActor_h +#define GateFluenceActor_h + +#include "G4Cache.hh" +#include "G4VPrimitiveScorer.hh" +#include "GateVActor.h" +#include "itkImage.h" +#include +#include + +namespace py = pybind11; + +class GateFluenceActor : public GateVActor { + +public: + // Constructor + GateFluenceActor(py::dict &user_info); + + virtual void ActorInitialize(); + + // Main function called every step in attached volume + virtual void SteppingAction(G4Step *); + + // Called every time a Run starts (all threads) + virtual void BeginOfRunAction(const G4Run *run); + + virtual void BeginOfEventAction(const G4Event *event); + + virtual void EndSimulationAction(); + + // Image type is 3D float by default + typedef itk::Image Image3DType; + + typedef itk::Image Image4DType; + typedef itk::Image ImageInt4DType; + using Size4DType = Image4DType::SizeType; + Size4DType size_4D; + + // The image is accessible on py side (shared by all threads) + Image3DType::Pointer cpp_fluence_image; + Image3DType::SizeType size_fluence; + double fVoxelVolume; + int NbOfThreads = 0; + + std::string fPhysicalVolumeName; + + G4ThreeVector fInitialTranslation; + std::string fHitType; +}; + +#endif // GateFluenceActor_h diff --git a/core/opengate_core/opengate_lib/GateGANSource.cpp b/core/opengate_core/opengate_lib/GateGANSource.cpp index 8e4e193b2..0b39f1fd8 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANSource.cpp @@ -175,6 +175,7 @@ void GateGANSource::GenerateOnePrimary(G4Event *event, // (if it is not set by the GAN, we may avoid to sample at each iteration) if (fPosition_is_set_by_GAN || fCurrentSkippedEvents == 0) position = GeneratePrimariesPosition(); + // FIXME change position is not set by GAN // direction // (if it is not set by the GAN, we may avoid to sample at each iteration) diff --git a/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp b/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp new file mode 100644 index 000000000..c21fe9e3e --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp @@ -0,0 +1,22 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "GateFluenceActor.h" + +void init_GateFluenceActor(py::module &m) { + py::class_, + GateVActor>(m, "GateFluenceActor") + .def(py::init()) + .def_readwrite("cpp_fluence_image", &GateFluenceActor::cpp_fluence_image) + .def_readwrite("fPhysicalVolumeName", + &GateFluenceActor::fPhysicalVolumeName); +} diff --git a/opengate/actors/actorbuilders.py b/opengate/actors/actorbuilders.py index d597ea94e..284d89a6e 100644 --- a/opengate/actors/actorbuilders.py +++ b/opengate/actors/actorbuilders.py @@ -1,5 +1,5 @@ from .arfactors import ARFActor, ARFTrainingDatasetActor -from .doseactors import DoseActor, LETActor +from .doseactors import DoseActor, LETActor, FluenceActor from .digitizers import ( DigitizerAdderActor, DigitizerReadoutActor, @@ -24,6 +24,7 @@ actor_type_names = { SimulationStatisticsActor, DoseActor, + FluenceActor, LETActor, SourceInfoActor, PhaseSpaceActor, diff --git a/opengate/actors/doseactors.py b/opengate/actors/doseactors.py index 292f63c7b..7c83a73d1 100644 --- a/opengate/actors/doseactors.py +++ b/opengate/actors/doseactors.py @@ -200,7 +200,6 @@ def StartSimulationAction(self): self.output_origin = self.user_info.output_origin def EndSimulationAction(self): - # print(lol) g4.GateDoseActor.EndSimulationAction(self) # Get the itk image from the cpp side @@ -515,3 +514,113 @@ def EndSimulationAction(self): itk.imwrite(self.py_numerator_image, ensure_filename_is_str(fPath)) fPath = fPath.replace("_numerator", "_denominator") itk.imwrite(self.py_denominator_image, ensure_filename_is_str(fPath)) + + +class FluenceActor(g4.GateFluenceActor, ActorBase): + """ + FluenceActor: compute a 3D map of fluence + + FIXME: add scatter order and uncertainty + """ + + type_name = "FluenceActor" + + def set_default_user_info(user_info): + ActorBase.set_default_user_info(user_info) + # required user info, default values + mm = g4_units.mm + user_info.size = [10, 10, 10] + user_info.spacing = [1 * mm, 1 * mm, 1 * mm] + user_info.output = "fluence.mhd" + user_info.translation = [0, 0, 0] + user_info.physical_volume_index = None + user_info.uncertainty = False + user_info.scatter = False + + def __init__(self, user_info): + ActorBase.__init__(self, user_info) + g4.GateFluenceActor.__init__(self, user_info.__dict__) + # attached physical volume (at init) + self.g4_phys_vol = None + # default image (py side) + self.py_fluence_image = None + + def __str__(self): + u = self.user_info + s = f'FluenceActor "{u.name}": dim={u.size} spacing={u.spacing} {u.output} tr={u.translation}' + return s + + def __getstate__(self): + # superclass getstate + DoseActor.__getstate__(self) + return self.__dict__ + + def initialize(self, volume_engine=None): + super().initialize(volume_engine) + # create itk image (py side) + size = np.array(self.user_info.size) + spacing = np.array(self.user_info.spacing) + self.py_fluence_image = create_3d_image(size, spacing) + # compute the center, using translation and half pixel spacing + self.img_origin_during_run = ( + -size * spacing / 2.0 + spacing / 2.0 + self.user_info.translation + ) + # for initialization during the first run + self.first_run = True + + def StartSimulationAction(self): + # init the origin and direction according to the physical volume + # (will be updated in the BeginOfRun) + attached_to_volume = self.volume_engine.get_volume(self.user_info.mother) + if self.user_info.physical_volume_index is None: + physical_volume_index = 0 + else: + physical_volume_index = self.user_info.physical_volume_index + try: + self.g4_phys_vol = attached_to_volume.g4_physical_volumes[ + physical_volume_index + ] + except IndexError: + fatal( + f"Error in the FluenceActor {self.user_info.name}. " + f"Could not find the physical volume with index {physical_volume_index} " + f"in volume '{self.user_info.mother}' to which this actor is attached. " + ) + align_image_with_physical_volume( + attached_to_volume, + self.py_fluence_image, + initial_translation=self.user_info.translation, + ) + + # Set the real physical volume name + self.fPhysicalVolumeName = str(self.g4_phys_vol.GetName()) + + # FIXME for multiple run and motion + if not self.first_run: + warning(f"Not implemented yet: FluenceActor with several runs") + # send itk image to cpp side, copy data only the first run. + update_image_py_to_cpp( + self.py_fluence_image, self.cpp_fluence_image, self.first_run + ) + + # now, indicate the next run will not be the first + self.first_run = False + + def EndSimulationAction(self): + g4.GateFluenceActor.EndSimulationAction(self) + + # Get the itk image from the cpp side + # Currently a copy. Maybe later as_pyarray ? + self.py_fluence_image = get_cpp_image(self.cpp_fluence_image) + + # set the property of the output image: + origin = self.img_origin_during_run + self.py_fluence_image.SetOrigin(origin) + + # write the image at the end of the run + # FIXME : maybe different for several runs + if self.user_info.output: + out_p = ensure_filename_is_str( + self.simulation.get_output_path(self.user_info.output) + ) + itk.imwrite(self.py_fluence_image, out_p) diff --git a/opengate/tests/data b/opengate/tests/data index 19376c4fa..4347b2502 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit 19376c4fa97c41c20f4b96b1e946cff0a1e5bfee +Subproject commit 4347b25021125285a7edbd9ad7d555be9f23e93c From 217cf0a932dfcd62feb534d193487d3c95173eb6 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 17:57:53 +0100 Subject: [PATCH 21/37] no options yet --- opengate/actors/doseactors.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/opengate/actors/doseactors.py b/opengate/actors/doseactors.py index 7c83a73d1..e021b8b31 100644 --- a/opengate/actors/doseactors.py +++ b/opengate/actors/doseactors.py @@ -567,6 +567,9 @@ def initialize(self, volume_engine=None): ) # for initialization during the first run self.first_run = True + # no options yet + if self.user_info.uncertainty or self.user_info.scatter: + fatal(f"FluenceActor : uncertainty and scatter not implemented yet") def StartSimulationAction(self): # init the origin and direction according to the physical volume From ee306a7c9466b5ce475ff145b4e218a1ec2df1f1 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 18:13:29 +0100 Subject: [PATCH 22/37] remove debug print --- opengate/engines.py | 7 ++++--- opengate/geometry/utility.py | 1 - opengate/managers.py | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/opengate/engines.py b/opengate/engines.py index 472f2853a..c68491310 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -1480,9 +1480,10 @@ def initialize_g4_verbose(self): else: # no Geant4 output ui = UIsessionSilent() - self.simulation.add_g4_command_after_init( - f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" - ) + if self.simulation.g4_verbose_level_tracking >= 0: + self.simulation.add_g4_command_after_init( + f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" + ) # it is also possible to set ui=None for 'default' output # we must keep a ref to ui_session self.ui_session = ui diff --git a/opengate/geometry/utility.py b/opengate/geometry/utility.py index 24b70fed4..3bec06d18 100644 --- a/opengate/geometry/utility.py +++ b/opengate/geometry/utility.py @@ -188,7 +188,6 @@ def get_transform_world_to_local(volume): ctr = volume.translation_list[i] crot = volume.rotation_list[i] for vol in volume.ancestor_volumes[::-1]: - print(vol.name) crot = np.matmul(vol.rotation_list[0], crot) ctr = vol.rotation_list[0].dot(ctr) + vol.translation_list[0] cumulative_translation.append(ctr) diff --git a/opengate/managers.py b/opengate/managers.py index 4231786d9..3a15c5afc 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -909,7 +909,7 @@ class Simulation(GateObject): ), "g4_verbose": (False, {"doc": "Switch on/off Geant4's verbose output."}), "g4_verbose_level_tracking": ( - 0, + -1, { "doc": "Activate verbose tracking in Geant4 via G4 command '/tracking/verbose g4_verbose_level_tracking'." }, From 418f32a84ced1ae064ba511216df76bfe1631aab Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 19:04:14 +0100 Subject: [PATCH 23/37] add data --- opengate/tests/data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/tests/data b/opengate/tests/data index 4347b2502..c7b2d4733 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit 4347b25021125285a7edbd9ad7d555be9f23e93c +Subproject commit c7b2d4733d018b1260c1b7e80c87d7c829e47c4d From b7214be226ca22454623a7a75149c64ca848fcbd Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 19:10:59 +0100 Subject: [PATCH 24/37] remove scatter option --- .../src/test067_cbct_fluence_actor_mt.py | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100755 opengate/tests/src/test067_cbct_fluence_actor_mt.py diff --git a/opengate/tests/src/test067_cbct_fluence_actor_mt.py b/opengate/tests/src/test067_cbct_fluence_actor_mt.py new file mode 100755 index 000000000..0d5b5bfa4 --- /dev/null +++ b/opengate/tests/src/test067_cbct_fluence_actor_mt.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "gate_test067", "test067") + + sim = gate.Simulation() + + ui = sim.user_info + # ui.visu = True + ui.visu_type = "vrml" + ui.random_engine = "MersenneTwister" + ui.random_seed = "auto" + ui.number_of_threads = 2 + + # units + m = gate.g4_units.m + cm = gate.g4_units.cm + nm = gate.g4_units.nm + cm3 = gate.g4_units.cm3 + keV = gate.g4_units.keV + MeV = gate.g4_units.MeV + mm = gate.g4_units.mm + Bq = gate.g4_units.Bq + gcm3 = gate.g4_units.g / cm3 + + # world + world = sim.world + world.size = [3 * m, 3 * m, 3 * m] + world.material = "G4_AIR" + + # waterbox + waterbox = sim.add_volume("Box", "waterbox") + waterbox.size = [20 * cm, 20 * cm, 20 * cm] + waterbox.translation = [0 * cm, 0 * cm, 15 * cm] + waterbox.material = "G4_WATER" + + # CBCT gantry source + gantry = sim.add_volume("Box", "CBCT_gantry") + gantry.size = [0.2 * m, 0.2 * m, 0.2 * m] + gantry.translation = [1000 * mm, 0, 0] + gantry.material = "G4_AIR" + gantry.color = [0, 1, 1, 1] + + # CBCT detector plane + detector_plane = sim.add_volume("Box", "CBCT_detector_plane") + detector_plane.size = [10 * mm, 1409.6 * mm, 1409.6 * mm] + detector_plane.translation = [-536 * mm, 0, 0] + detector_plane.material = "G4_AIR" + detector_plane.color = [1, 0, 0, 1] + + # physics + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option1" + sim.physics_manager.set_production_cut("world", "all", 10 * mm) + + # actor + detector_actor = sim.add_actor("FluenceActor", "detector_actor") + detector_actor.mother = detector_plane.name + detector_actor.output = paths.output / "fluence.mhd" + detector_actor.spacing = [10 * mm, 5 * mm, 5 * mm] + detector_actor.size = [1, 100, 100] + + # source + source = sim.add_source("GenericSource", "mysource") + source.mother = gantry.name + source.particle = "gamma" + source.energy.mono = 60 * keV + source.position.type = "box" + source.position.size = [1 * nm, 2 * 8 * mm, 2 * 8 * mm] + source.direction.type = "focused" + # FIXME warning in world coord ! Should be in mother coord system + source.direction.focus_point = [gantry.translation[0] - 60 * mm, 0, 0] + source.n = 50000 / ui.number_of_threads + + if ui.visu: + source.n = 100 + + # statistics + stats = sim.add_actor("SimulationStatisticsActor", "stats") + stats.track_types_flag = True + stats.output = paths.output / "stats.txt" + + # run + output = sim.run() + + # print output statistics + stats = output.get_actor("stats") + print(stats) + + # check images + is_ok = utility.assert_images( + paths.gate_output / "detector.mhd", + paths.output / "fluence.mhd", + stats, + tolerance=44, + axis="y", + ) + + utility.test_ok(is_ok) From 57ecd5c5a5a17f0325f27d7d9ed1402aaac59a77 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 4 Dec 2023 21:46:09 +0100 Subject: [PATCH 25/37] arf plane side --- core/opengate_core/opengate_lib/GateARFActor.cpp | 5 ++++- core/opengate_core/opengate_lib/GateARFActor.h | 1 + opengate/actors/arfactors.py | 7 ++++--- opengate/tests/src/test043_garf.py | 7 ++++--- opengate/tests/src/test043_garf_analog.py | 3 +-- opengate/tests/src/test043_garf_mt.py | 1 + .../tests/src/test043_garf_training_dataset_full_wip.py | 3 +++ opengate/tests/src/test047_gan_vox_source_cond.py | 9 +++++---- 8 files changed, 23 insertions(+), 13 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateARFActor.cpp b/core/opengate_core/opengate_lib/GateARFActor.cpp index ad9ec34ee..964924175 100644 --- a/core/opengate_core/opengate_lib/GateARFActor.cpp +++ b/core/opengate_core/opengate_lib/GateARFActor.cpp @@ -16,6 +16,7 @@ GateARFActor::GateARFActor(py::dict &user_info) : GateVActor(user_info, true) { fActions.insert("EndOfRunAction"); // User option: batch size fBatchSize = DictGetInt(user_info, "batch_size"); + fKeepNegativeSide = DictGetBool(user_info, "flip_plane"); } void GateARFActor::SetARFFunction(ARFFunctionType &f) { fApply = f; } @@ -50,7 +51,9 @@ void GateARFActor::SteppingAction(G4Step *step) { dir = pre->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(dir); // which side of the plane ? - if (dir[2] < 0) + if (!fKeepNegativeSide && dir[2] < 0) + return; + if (fKeepNegativeSide && dir[2] > 0) return; l.fCurrentNumberOfHits++; diff --git a/core/opengate_core/opengate_lib/GateARFActor.h b/core/opengate_core/opengate_lib/GateARFActor.h index 7f0902858..3ff7f66f1 100644 --- a/core/opengate_core/opengate_lib/GateARFActor.h +++ b/core/opengate_core/opengate_lib/GateARFActor.h @@ -54,6 +54,7 @@ class GateARFActor : public GateVActor { protected: int fBatchSize; ARFFunctionType fApply; + bool fKeepNegativeSide; // For MT, all threads local variables are gathered here struct threadLocalT { diff --git a/opengate/actors/arfactors.py b/opengate/actors/arfactors.py index 50fd340bd..01e796244 100644 --- a/opengate/actors/arfactors.py +++ b/opengate/actors/arfactors.py @@ -105,6 +105,7 @@ def set_default_user_info(user_info): user_info.verbose_batch = False user_info.output = "" user_info.enable_hit_slice = False + user_info.flip_plane = False # Can be cpu / auto / gpu user_info.gpu_mode = "auto" @@ -300,6 +301,6 @@ def EndSimulationAction(self): self.output_image, ensure_filename_is_str(self.user_info.output) ) - # FIXME debug - print(f"{self.debug_nb_hits_before=}") - print(f"{self.debug_nb_hits=}") + # debug + # print(f"{self.debug_nb_hits_before=}") + # print(f"{self.debug_nb_hits=}") diff --git a/opengate/tests/src/test043_garf.py b/opengate/tests/src/test043_garf.py index 2fad11a53..0a3b3fb5a 100755 --- a/opengate/tests/src/test043_garf.py +++ b/opengate/tests/src/test043_garf.py @@ -69,6 +69,8 @@ arf.distance_to_crystal = crystal_dist # 74.625 * mm arf.distance_to_crystal = 74.625 * mm ## FIXME arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" + # arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v034.pth" + arf.flip_plane = True # because the training was backside arf.enable_hit_slice = True # the gpu_mode should be "auto" but when running on github actions, # with set "cpu" to avoid mps errors @@ -77,8 +79,7 @@ s = sim.add_actor("SimulationStatisticsActor", "stats") s.track_types_flag = True - # FIXME - """# start simulation to check if ok with start_new_process + # start simulation to check if ok with start_new_process s1a = s1.activity s2a = s2.activity s3a = s3.activity @@ -90,7 +91,7 @@ # restart simulation s1.activity = s1a s2.activity = s2a - s3.activity = s3a""" + s3.activity = s3a sim.run(False) # print results at the end diff --git a/opengate/tests/src/test043_garf_analog.py b/opengate/tests/src/test043_garf_analog.py index 8ea3a5b63..8e7671bad 100755 --- a/opengate/tests/src/test043_garf_analog.py +++ b/opengate/tests/src/test043_garf_analog.py @@ -36,8 +36,7 @@ spect, crystal = gate_spect.add_ge_nm67_spect_head( sim, "spect", collimator_type=colli, debug=sim.visu ) - spect_translation = 15 * cm - spect.translation = [0, 0, -spect_translation] + spect.translation = [0, 0, -15 * cm] crystal_name = f"{spect.name}_crystal" # physics diff --git a/opengate/tests/src/test043_garf_mt.py b/opengate/tests/src/test043_garf_mt.py index fed7555d6..04bba5746 100755 --- a/opengate/tests/src/test043_garf_mt.py +++ b/opengate/tests/src/test043_garf_mt.py @@ -67,6 +67,7 @@ arf.distance_to_crystal = 74.625 * mm arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" arf.enable_hit_slice = True + arf.flip_plane = True # because the training was backside arf.gpu_mode = ( utility.get_gpu_mode_for_tests() ) # should be "auto" but "cpu" for macOS github actions to avoid mps errors diff --git a/opengate/tests/src/test043_garf_training_dataset_full_wip.py b/opengate/tests/src/test043_garf_training_dataset_full_wip.py index 8c6b446fe..f8de845f8 100755 --- a/opengate/tests/src/test043_garf_training_dataset_full_wip.py +++ b/opengate/tests/src/test043_garf_training_dataset_full_wip.py @@ -96,3 +96,6 @@ print(stat) skip = gate.sources.generic.get_source_skipped_events(sim.output, "s1") print(f"Nb of skip particles {skip} {(skip / stat.counts.event_count) * 100:.2f}%") + + # garf + # garf_train train_arf_v034.json /home/dsarrut/src/gate2/opengate/opengate/tests/src/../output/test043_arf_training_dataset_large.root a.pth diff --git a/opengate/tests/src/test047_gan_vox_source_cond.py b/opengate/tests/src/test047_gan_vox_source_cond.py index c26c978ab..14072fd16 100755 --- a/opengate/tests/src/test047_gan_vox_source_cond.py +++ b/opengate/tests/src/test047_gan_vox_source_cond.py @@ -20,7 +20,8 @@ # main options sim.g4_verbose = False sim.g4_verbose_level = 1 - sim.visu = False + # sim.visu = True + sim.visu_type = "vrml" sim.number_of_threads = 1 sim.random_seed = 123456789 activity_bq = 1e6 @@ -77,9 +78,9 @@ source = sim.add_source("GANSource", "source") source.mother = "ct" source.cond_image = paths.data / "source_three_areas_crop_3.5mm.mhd" - source.position.translation = gate.image.get_translation_between_images_center( - str(ct.image), str(source.cond_image) - ) + # source.position.translation = gate.image.get_translation_between_images_center( + # str(ct.image), str(source.cond_image) + # ) source.particle = "alpha" source.activity = activity_bq * Bq / sim.number_of_threads source.compute_directions = True From dffd33ca85a45e6b164b60d36b321bc44afcffdd Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 5 Dec 2023 09:24:04 +0100 Subject: [PATCH 26/37] update gaga translation --- opengate/tests/src/test043_garf.py | 4 ++-- opengate/tests/src/test043_garf_mt.py | 3 ++- opengate/tests/src/test047_gan_vox_source_cond.py | 8 +++++--- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/opengate/tests/src/test043_garf.py b/opengate/tests/src/test043_garf.py index 0a3b3fb5a..5dc16581c 100755 --- a/opengate/tests/src/test043_garf.py +++ b/opengate/tests/src/test043_garf.py @@ -68,8 +68,8 @@ arf.verbose_batch = True arf.distance_to_crystal = crystal_dist # 74.625 * mm arf.distance_to_crystal = 74.625 * mm ## FIXME - arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" - # arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v034.pth" + # arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" + arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v034.pth" arf.flip_plane = True # because the training was backside arf.enable_hit_slice = True # the gpu_mode should be "auto" but when running on github actions, diff --git a/opengate/tests/src/test043_garf_mt.py b/opengate/tests/src/test043_garf_mt.py index 04bba5746..025108c4c 100755 --- a/opengate/tests/src/test043_garf_mt.py +++ b/opengate/tests/src/test043_garf_mt.py @@ -65,7 +65,8 @@ arf.verbose_batch = True arf.distance_to_crystal = crystal_dist # 74.625 * mm arf.distance_to_crystal = 74.625 * mm - arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" + # arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v3.pth" + arf.pth_filename = test43.paths.gate_data / "pth" / "arf_Tc99m_v034.pth" arf.enable_hit_slice = True arf.flip_plane = True # because the training was backside arf.gpu_mode = ( diff --git a/opengate/tests/src/test047_gan_vox_source_cond.py b/opengate/tests/src/test047_gan_vox_source_cond.py index 14072fd16..ade534276 100755 --- a/opengate/tests/src/test047_gan_vox_source_cond.py +++ b/opengate/tests/src/test047_gan_vox_source_cond.py @@ -78,9 +78,11 @@ source = sim.add_source("GANSource", "source") source.mother = "ct" source.cond_image = paths.data / "source_three_areas_crop_3.5mm.mhd" - # source.position.translation = gate.image.get_translation_between_images_center( - # str(ct.image), str(source.cond_image) - # ) + source.position.translation = gate.image.get_translation_between_images_center( + str(ct.image), str(source.cond_image) + ) + source.position.translation = source.position.translation / 2.0 + print(f"translation {source.position.translation}") source.particle = "alpha" source.activity = activity_bq * Bq / sim.number_of_threads source.compute_directions = True From b54b86b308788d5496322d2a34c265334c49fa59 Mon Sep 17 00:00:00 2001 From: David Date: Wed, 6 Dec 2023 14:14:42 +0100 Subject: [PATCH 27/37] typo --- opengate/exception.py | 1 - opengate/sources/gansources.py | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/opengate/exception.py b/opengate/exception.py index 9d90c5ada..7c4b1db2f 100644 --- a/opengate/exception.py +++ b/opengate/exception.py @@ -1,6 +1,5 @@ import inspect import colored -import sys import opengate_core as g4 from .logger import log diff --git a/opengate/sources/gansources.py b/opengate/sources/gansources.py index 63d3dcfe3..c7231688f 100644 --- a/opengate/sources/gansources.py +++ b/opengate/sources/gansources.py @@ -577,7 +577,7 @@ def generator(self, source): print(f"Generate {n} particles from GAN ", end="") # generate samples (this is the most time-consuming part) - fake = self.gaga.generate_samples_no_cond( + fake = self.gaga.generate_samples_non_cond( g.params, g.G, n=n, @@ -747,7 +747,7 @@ def generator(self, source): print(f"Generate {n} particles from GAN ", end="") # generate samples (this is the most time-consuming part) - fake = self.gaga.generate_samples_no_cond( + fake = self.gaga.generate_samples_non_cond( g.params, g.G, n=n, From cfea8e67f10df2447d22524d7bd6135e027944c8 Mon Sep 17 00:00:00 2001 From: David Date: Wed, 6 Dec 2023 14:41:52 +0100 Subject: [PATCH 28/37] no nb thread --- opengate/actors/miscactors.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/opengate/actors/miscactors.py b/opengate/actors/miscactors.py index ead2f3c44..5b12f1566 100644 --- a/opengate/actors/miscactors.py +++ b/opengate/actors/miscactors.py @@ -71,14 +71,6 @@ def sps(self): return self.counts.step_count / self.counts.duration * sec return 0 - @property - def nb_thread(self): - if self.simulation is not None: - thread = self.simulation.number_of_threads - else: - thread = "?" - return thread - @property def simu_start_time(self): if not self.simulation is None: From 2e1b4a6c62d6f234092fcd3a10db48f9ad4e1fb7 Mon Sep 17 00:00:00 2001 From: David Date: Mon, 11 Dec 2023 15:13:33 +0100 Subject: [PATCH 29/37] doserate with visu --- opengate/actors/doseactors.py | 2 +- opengate/bin/dose_rate | 4 +- opengate/contrib/dose/doserate.py | 51 ++++++++++++++----- opengate/engines.py | 7 +-- .../src/test067_cbct_fluence_actor_mt.py | 4 +- 5 files changed, 47 insertions(+), 21 deletions(-) diff --git a/opengate/actors/doseactors.py b/opengate/actors/doseactors.py index e021b8b31..947428b2c 100644 --- a/opengate/actors/doseactors.py +++ b/opengate/actors/doseactors.py @@ -508,7 +508,7 @@ def EndSimulationAction(self): ) itk.imwrite(self.py_LETd_image, ensure_filename_is_str(fPath)) - # for parrallel computation we need to provide both outputs + # for parallel computation we need to provide both outputs if self.user_info.separate_output: fPath = fPath.replace(".mhd", "_numerator.mhd") itk.imwrite(self.py_numerator_image, ensure_filename_is_str(fPath)) diff --git a/opengate/bin/dose_rate b/opengate/bin/dose_rate index 90327acd0..fcfe99785 100755 --- a/opengate/bin/dose_rate +++ b/opengate/bin/dose_rate @@ -35,10 +35,10 @@ def go(json_param, output_folder): sim = create_simulation(param) # run - output = sim.start() + sim.run() # print results at the end - stats = output.get_actor("Stats") + stats = sim.output.get_actor("Stats") print(stats) stats.write(param.output_folder / "stats.txt") print(f"Output in {param.output_folder}") diff --git a/opengate/contrib/dose/doserate.py b/opengate/contrib/dose/doserate.py index 6b78a278d..7307fc0a1 100755 --- a/opengate/contrib/dose/doserate.py +++ b/opengate/contrib/dose/doserate.py @@ -10,6 +10,19 @@ def create_simulation(param): + """ + param is dict with: + - output_folder + - visu + - number_of_threads + - ct_image + - density_tolerance_gcm3 + - table_mat + - table_density + - verbose: + - radionuclide + - activity_bq + """ # create the simulation sim = Simulation() @@ -17,6 +30,7 @@ def create_simulation(param): ui = sim.user_info ui.g4_verbose = False ui.visu = param.visu + ui.visu_type = "vrml" ui.number_of_threads = param.number_of_threads ui.verbose_level = INFO @@ -34,17 +48,27 @@ def create_simulation(param): world.size = [2 * m, 2 * m, 2 * m] # CT image - ct = sim.add_volume("Image", "ct") - ct.image = param.ct_image - ct.material = "G4_AIR" # material used by default - tol = param.density_tolerance_gcm3 * gcm3 - ct.voxel_materials, materials = HounsfieldUnit_to_material( - sim, tol, param.table_mat, param.table_density - ) - if param.verbose: - print(f'Density tolerance = {g4_best_unit(tol, "Volumic Mass")}') - print(f"Number of materials in the CT : {len(ct.voxel_materials)} materials") - ct.dump_label_image = param.output_folder / "labels.mhd" + if ui.visu: + ct = sim.add_volume("Box", "ct") + info = read_image_info(param.ct_image) + ct.size = info.size + ct.material = "G4_WATER" + ct.color = [0, 0, 1, 1] + ui.number_of_threads = 1 + else: + ct = sim.add_volume("Image", "ct") + ct.image = param.ct_image + ct.material = "G4_AIR" # material used by default + tol = param.density_tolerance_gcm3 * gcm3 + ct.voxel_materials, materials = HounsfieldUnit_to_material( + sim, tol, param.table_mat, param.table_density + ) + if param.verbose: + print(f'Density tolerance = {g4_best_unit(tol, "Volumic Mass")}') + print( + f"Number of materials in the CT : {len(ct.voxel_materials)} materials" + ) + ct.dump_label_image = param.output_folder / "labels.mhd" # some radionuclides choice # (user of this function can still change @@ -74,7 +98,7 @@ def create_simulation(param): # cuts sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4" - sim.physics_manager.enable_decay = True # FIXME + sim.physics_manager.enable_decay = True sim.physics_manager.set_production_cut("world", "all", 1 * m) sim.physics_manager.set_production_cut("ct", "all", 1 * mm) @@ -88,7 +112,8 @@ def create_simulation(param): # translate the dose the same way as the source dose.translation = source.position.translation # set the origin of the dose like the source - dose.img_coord_system = True + if not ui.visu: + dose.img_coord_system = True dose.hit_type = "random" dose.uncertainty = False diff --git a/opengate/engines.py b/opengate/engines.py index e6f591f08..4a03947c7 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -1245,9 +1245,10 @@ def initialize_g4_verbose(self): else: # no Geant4 output ui = UIsessionSilent() - self.simulation.add_g4_command_after_init( - f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" - ) + if self.simulation.g4_verbose_level_tracking != -1: + self.simulation.add_g4_command_after_init( + f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" + ) # it is also possible to set ui=None for 'default' output # we must keep a ref to ui_session self.ui_session = ui diff --git a/opengate/tests/src/test067_cbct_fluence_actor_mt.py b/opengate/tests/src/test067_cbct_fluence_actor_mt.py index 0d5b5bfa4..58f903a9e 100755 --- a/opengate/tests/src/test067_cbct_fluence_actor_mt.py +++ b/opengate/tests/src/test067_cbct_fluence_actor_mt.py @@ -84,10 +84,10 @@ stats.output = paths.output / "stats.txt" # run - output = sim.run() + sim.run() # print output statistics - stats = output.get_actor("stats") + stats = sim.output.get_actor("stats") print(stats) # check images From 4586847c2432919ed7140ccb9a98445dcab5e9d2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 25 Dec 2023 20:35:50 +0000 Subject: [PATCH 30/37] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/psf/black: 23.12.0 → 23.12.1](https://github.com/psf/black/compare/23.12.0...23.12.1) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6abf66756..0165d07d4 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -5,7 +5,7 @@ repos: - id: end-of-file-fixer - id: trailing-whitespace - repo: https://github.com/psf/black - rev: 23.12.0 + rev: 23.12.1 hooks: - id: black - repo: https://github.com/pre-commit/mirrors-clang-format From 8c73b5738721a8aabbc9004ec69277f27782d6a0 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 8 Jan 2024 11:59:29 +0100 Subject: [PATCH 31/37] Change minimal version of gaga-phsp and garf Ignore tests test066_spect_gaga_garf_4_analyse1.py and test066_spect_gaga_garf_5_analyse2.py because for the moment they lead to errors. They need output data from previous tests and they do not create data, just compare them --- .github/workflows/main.yml | 3 +-- opengate/actors/arfactors.py | 2 +- opengate/bin/opengate_tests | 2 ++ opengate/sources/gansources.py | 2 +- 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index bff6f8a1f..e9e05a1b1 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -363,8 +363,7 @@ jobs: pip install torch fi pip install SimpleITK - pip install gaga_phsp - pip install garf + pip install gaga_phsp==0.7.1 pip install dist/opengate_core-*-${PYTHONFOLDER}-${OSNAME}*_${PLATFORM}64.whl pip install dist/opengate-*.whl export GIT_SSL_NO_VERIFY=1 diff --git a/opengate/actors/arfactors.py b/opengate/actors/arfactors.py index 01e796244..ef7dec418 100644 --- a/opengate/actors/arfactors.py +++ b/opengate/actors/arfactors.py @@ -30,7 +30,7 @@ def import_garf(): from packaging import version garf_version = pkg_resources.get_distribution("garf").version - garf_minimal_version = "2.4" + garf_minimal_version = "2.5" if version.parse(garf_version) < version.parse(garf_minimal_version): fatal( "The minimal version of garf is not correct. You should install at least the version " diff --git a/opengate/bin/opengate_tests b/opengate/bin/opengate_tests index 599f3a0c5..85364624f 100755 --- a/opengate/bin/opengate_tests +++ b/opengate/bin/opengate_tests @@ -64,6 +64,8 @@ def go(test_id, random_tests): "test066_spect_gaga_garf_1_reference.py", # ignored because reference data (too long) "test066_spect_gaga_garf_2.py", # ignored because reference data (too long, GPU) "test066_spect_gaga_garf_3_standalone.py", # ignored because too long (GPU) + "test066_spect_gaga_garf_4_analyse1.py", + "test066_spect_gaga_garf_5_analyse2.py", ] onlyfiles = [ diff --git a/opengate/sources/gansources.py b/opengate/sources/gansources.py index c7231688f..5bae26db5 100644 --- a/opengate/sources/gansources.py +++ b/opengate/sources/gansources.py @@ -35,7 +35,7 @@ def import_gaga_phsp(): from packaging import version gaga_version = pkg_resources.get_distribution("gaga_phsp").version - gaga_minimal_version = "0.7.0" + gaga_minimal_version = "0.7.1" if version.parse(gaga_version) < version.parse(gaga_minimal_version): fatal( "The minimal version of gaga_phsp is not correct. You should install at least the version " From 4b7055cd83079dd22229d8cfb2fa15279d66a73c Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 8 Jan 2024 14:55:22 +0100 Subject: [PATCH 32/37] Update data --- opengate/tests/data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/tests/data b/opengate/tests/data index a0e299e40..3d8000895 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit a0e299e40bb878b7bdeed76d0f9f1a6bfa993557 +Subproject commit 3d800089595f329ab46cedc45554729e5df43c44 From 5aa323e549e51c08e1a7aa5b888e2d350ab3d9c1 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Mon, 8 Jan 2024 14:55:46 +0100 Subject: [PATCH 33/37] Try newest version of macos fir github actions --- .github/workflows/main.yml | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index e9e05a1b1..d61d7fb65 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -17,10 +17,10 @@ jobs: strategy: fail-fast: false matrix: - os: [ubuntu-latest, macos-11, windows-latest] + os: [ubuntu-latest, macos-latest, windows-latest] python-version: [3.8, 3.9, '3.10', '3.11'] exclude: - - os: macos-11 + - os: macos-latest python-version: '3.11' env: @@ -36,7 +36,7 @@ jobs: if [[ ${{ matrix.os }} == "ubuntu-latest" ]]; then export GIT_SSL_NO_VERIFY=1 git submodule update --init --recursive - elif [[ ${{ matrix.os }} == "macos-11" ]]; then + elif [[ ${{ matrix.os }} == "macos-latest" ]]; then export GIT_SSL_NO_VERIFY=1 git submodule update --init --recursive else @@ -56,7 +56,7 @@ jobs: varOS=`cat /etc/os-release | grep "VERSION=" | grep -oP '(?<=\").*?(?=\")'` varOS=($varOS) echo "release=${varOS[0]}" >> $GITHUB_OUTPUT - elif [[ ${{ matrix.os }} == "macos-11" ]]; then + elif [[ ${{ matrix.os }} == "macos-latest" ]]; then varOS=`sw_vers | grep "ProductVersion:"` varOS="${varOS#*:}" echo "release=${varOS:1}" >> $GITHUB_OUTPUT @@ -106,22 +106,22 @@ jobs: mv dist_opengate/* dist/ fi - uses: conda-incubator/setup-miniconda@v2 - if: (matrix.os == 'macos-11') || (matrix.os == 'windows-latest') + if: (matrix.os == 'macos-latest') || (matrix.os == 'windows-latest') with: auto-update-conda: true activate-environment: opengate_core python-version: ${{ matrix.python-version }} - name: Set up Homebrew - if: matrix.os == 'macos-11' + if: matrix.os == 'macos-latest' id: set-up-homebrew uses: Homebrew/actions/setup-homebrew@master - name: Create opengate_core Wheel Mac - if: matrix.os == 'macos-11' + if: matrix.os == 'macos-latest' shell: bash -l {0} run: | brew update - rm -rf /usr/local/bin/python3.1*-config /usr/local/bin/2to3-3.1* /usr/local/bin/idle3.1* /usr/local/bin/pydoc3.1* /usr/local/bin/python3.1* - rm -rf /usr/local/bin/python3-config /usr/local/bin/2to3 /usr/local/bin/idle3 /usr/local/bin/pydoc3 /usr/local/bin/python3 + #rm -rf /usr/local/bin/python3.1*-config /usr/local/bin/2to3-3.1* /usr/local/bin/idle3.1* /usr/local/bin/pydoc3.1* /usr/local/bin/python3.1* + #rm -rf /usr/local/bin/python3-config /usr/local/bin/2to3 /usr/local/bin/idle3 /usr/local/bin/pydoc3 /usr/local/bin/python3 brew install --force --verbose --overwrite \ ccache \ fftw \ From b99e6001ceb2847ec5026ca0c8f358a70d1df8f2 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Tue, 9 Jan 2024 17:05:11 +0100 Subject: [PATCH 34/37] Support minimal version of python 3.8 --- core/setup.py | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/core/setup.py b/core/setup.py index fbe53e276..981d9f89e 100755 --- a/core/setup.py +++ b/core/setup.py @@ -177,7 +177,7 @@ def build_extension(self, ext): packages=find_packages(), package_data=package_data, zip_safe=False, - python_requires=">=3.5", + python_requires=">=3.8", include_package_data=True, classifiers=( "Programming Language :: Python :: 3", diff --git a/setup.py b/setup.py index a78e850bf..8615e7fe1 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ long_description_content_type="text/markdown", url="https://github.com/OpenGATE/opengate", packages=selected_packages, - python_requires=">=3.7", + python_requires=">=3.8", include_package_data=True, classifiers=( "Programming Language :: Python :: 3", From 098518b391475fe2e0576ec05df9cecb7ed66df6 Mon Sep 17 00:00:00 2001 From: tbaudier Date: Wed, 10 Jan 2024 09:30:03 +0100 Subject: [PATCH 35/37] Update cmake_prefix_path in the doc For ITK compilation, the build folder was not comprehensive --- docs/source/developer_guide.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/developer_guide.md b/docs/source/developer_guide.md index 18b37dee2..665c46658 100644 --- a/docs/source/developer_guide.md +++ b/docs/source/developer_guide.md @@ -58,8 +58,8 @@ For **ITK**, you need to compile with the following options: ```bash git clone --branch v5.2.1 https://github.com/InsightSoftwareConsortium/ITK.git --depth 1 -mkdir build-v5.2.1 -cd build-v5.2.1 +mkdir itk-build +cd itk-build cmake -DCMAKE_CXX_FLAGS=-std=c++17 \ -DBUILD_TESTING=OFF \ -DITK_USE_FFTWD=ON \ @@ -76,7 +76,7 @@ Once it is done, you can compile `opengate_core`. ```bash pip install colored cd /core -export CMAKE_PREFIX_PATH=/geant4.11-build/:/build-v5.1.0/:${CMAKE_PREFIX_PATH} +export CMAKE_PREFIX_PATH=/geant4.11-build/:/itk-build/:${CMAKE_PREFIX_PATH} pip install -e . -v ``` From a497cf4db04a79b426b328d32cfbed3106356d6d Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Thu, 11 Jan 2024 09:28:00 +0100 Subject: [PATCH 36/37] Change default output folder management See https://github.com/OpenGATE/opengate/issues/331 --- opengate/managers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opengate/managers.py b/opengate/managers.py index 9d32a1076..75c8fe7b0 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -1017,7 +1017,7 @@ class Simulation(GateObject): }, ), "output_dir": ( - "./output", + ".", { "doc": "Directory to which any output is written, " "unless an absolute path is provided for a specific output." @@ -1196,7 +1196,7 @@ def get_output_path(self, path=None, is_file_or_directory="file"): n = len(p_out.parts) - 1 # last item is the filename elif is_file_or_directory in ["dir", "Dir", "directory", "d"]: n = len(p_out.parts) # all items are part of the directory - if len(p_out.parts) > 0: + if len(p_out.parts) > 0 and n > 0: directory = Path(p_out.parts[0]) for i in range(n - 1): directory /= p_out.parts[i + 1] From 7aacedac4ae69a0c07eb2d9972633a89eabe5ea4 Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Thu, 11 Jan 2024 13:23:19 +0100 Subject: [PATCH 37/37] Towards version 10.0beta07 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index acea0ac76..e8a90296c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -10.0beta06 +10.0beta07