diff --git a/docs/source/user_guide/user_guide_how_to_convert_example_1.rst b/docs/source/user_guide/user_guide_how_to_convert_example_1.rst index 11356499f2..971a7a80a0 100644 --- a/docs/source/user_guide/user_guide_how_to_convert_example_1.rst +++ b/docs/source/user_guide/user_guide_how_to_convert_example_1.rst @@ -2,7 +2,7 @@ From GATE 9 to 10 - Example Proton CT ************************************* -This is an example of a proton beam with a spiral phantom placed in between two proton detectors. The spiral phantom was used in the article `Filtered backprojection proton CT reconstruction along most likely paths by Rit et al `_. The data generated by this simulation can be processed by the `PCT software `_. +This is an example of a proton beam with a spiral phantom placed in between two proton detectors. The spiral phantom was used in the article `Filtered backprojection proton CT reconstruction along most likely paths by Rit et al `_. A complete example of GATE proton CT simulation can be found `here `_, which can be run using the `protonct `_ executable. The data generated by this simulation can be processed by the `PCT software `_. .. image:: ../figures/proton_ct.png diff --git a/opengate/bin/protonct.py b/opengate/bin/protonct.py new file mode 100755 index 0000000000..1ac5a2eec1 --- /dev/null +++ b/opengate/bin/protonct.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import click +from opengate.contrib.protonct.protonct import protonct + +CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"]) + + +@click.command(context_settings=CONTEXT_SETTINGS) +@click.option("--output", "-o", default=None, help="output folder") +@click.option("--projections", "-p", default=720, help="number of projections") +@click.option( + "--protons-per-projection", + "-n", + default=1000, + help="number of protons per projection", +) +@click.option("--seed", "-s", type=int, help="random number generator seed") +@click.option("--visu", is_flag=True, help="Enable visualization") +@click.option("--verbose", is_flag=True, help="verbose execution") +def go(output, projections, protons_per_projection, seed, visu, verbose): + protonct(output, projections, protons_per_projection, seed, visu, verbose) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/opengate/contrib/protonct/__init__.py b/opengate/contrib/protonct/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/opengate/contrib/protonct/protonct.py b/opengate/contrib/protonct/protonct.py new file mode 100644 index 0000000000..e31e9904f0 --- /dev/null +++ b/opengate/contrib/protonct/protonct.py @@ -0,0 +1,166 @@ +import math +from scipy.spatial.transform import Rotation + +import opengate as gate + + +def protonct( + output, + projections=720, + protons_per_projection=1000, + seed=None, + visu=False, + verbose=False, +): + + # Units + nm = gate.g4_units.nm + mm = gate.g4_units.mm + cm = gate.g4_units.cm + m = gate.g4_units.m + sec = gate.g4_units.second + MeV = gate.g4_units.MeV + Bq = gate.g4_units.Bq + + # Simulation + sim = gate.Simulation() + + sim.random_engine = "MersenneTwister" + sim.random_seed = "auto" if seed is None else seed + sim.check_volumes_overlap = False + sim.visu = visu + sim.visu_type = "vrml" + sim.g4_verbose = False + sim.progress_bar = verbose + sim.number_of_threads = 1 + + sim.run_timing_intervals = [[j * sec, (j + 1) * sec] for j in range(projections)] + + # Misc + yellow = [1, 1, 0, 0.5] + blue = [0, 0, 1, 0.5] + + # Geometry + sim.volume_manager.add_material_database( + gate.utility.get_contrib_path() / "GateMaterials.db" + ) + sim.world.material = "Vacuum" + sim.world.size = [4 * m, 4 * m, 4 * m] + sim.world.color = [0, 0, 0, 0] + + # Phantom + + def add_spiral(sim): + # Mother of all + spiral = sim.add_volume("Tubs", name="Spiral") + spiral.rmin = 0 * cm + spiral.rmax = 10 * cm + spiral.dz = 40 * cm + spiral.material = "Water" + spiral.color = blue + spiral.rotation = Rotation.from_euler("yz", [90, 90], degrees=True).as_matrix() + + # Spiral rotation + tr, rot = gate.geometry.utility.volume_orbiting_transform( + "y", 0, 360, projections, spiral.translation, spiral.rotation + ) + spiral.add_dynamic_parametrisation(translation=tr, rotation=rot) + + # Spiral inserts + sradius = 4 + radius = list(range(0, 100 - sradius // 2, sradius)) + sangle = 139 + angles = [math.radians(a) for a in range(0, sangle * len(radius), sangle)] + posx = [radius[i] * math.cos(angles[i]) for i in range(len(radius))] + posy = [radius[i] * math.sin(angles[i]) for i in range(len(radius))] + + def add_spiral_insert( + sim, + mother, + name, + rmin=0 * mm, + rmax=1 * mm, + dz=40 * cm, + material="Aluminium", + translation=None, + color=None, + ): + + if translation is None: + translation = [0 * mm, 0 * mm, 0 * mm] + + if color is None: + color = yellow + + spiral_insert = sim.add_volume("Tubs", name=name) + spiral_insert.mother = mother.name + spiral_insert.rmin = rmin + spiral_insert.rmax = rmax + spiral_insert.dz = dz + spiral_insert.material = material + spiral_insert.translation = translation + spiral_insert.color = color + + for i in range(len(radius)): + add_spiral_insert( + sim, + spiral, + f"SpiralInsert{i:02d}", + translation=[posx[i] * mm, posy[i] * mm, 0], + ) + + add_spiral(sim) + + # Beam + source = sim.add_source("GenericSource", "mybeam") + source.particle = "proton" + source.energy.mono = 200 * MeV + source.energy.type = "mono" + source.position.type = "box" + source.position.size = [16 * mm, 1 * nm, 1 * nm] + source.position.translation = [0, 0, -1060 * mm] + source.direction.type = "focused" + source.direction.focus_point = [0, 0, -1000 * mm] + + if sim.visu: + # For visualisation speed, the number of particle is decreased + source.activity = 10 * Bq + else: + source.activity = protons_per_projection * Bq + + # Physics list + sim.physics_manager.physics_list_name = "QGSP_BIC_EMZ" + + # Phase spaces + + def add_detector(sim, name, translation): + plane = sim.add_volume("Box", "PlanePhaseSpace" + name) + plane.size = [400 * mm, 400 * mm, 1 * nm] + plane.translation = translation + plane.material = "Vacuum" + plane.color = yellow + + phase_space = sim.add_actor("PhaseSpaceActor", "PhaseSpace" + name) + phase_space.attached_to = plane.name + phase_space.attributes = [ + "RunID", + "EventID", + "TrackID", + "KineticEnergy", + "LocalTime", + "Position", + "Direction", + ] + phase_space.output_filename = f"{output}/PhaseSpace{name}.root" + + F = gate.actors.filters.GateFilterBuilder() + phase_space.filter = F.ParticleName == "proton" + + add_detector(sim, "In", [0, 0, -110 * mm]) + add_detector(sim, "Out", [0, 0, 110 * mm]) + + # Particle stats + stat = sim.add_actor("SimulationStatisticsActor", "stat") + stat.output_filename = f"{output}/protonct.txt" + + sim.run() diff --git a/opengate/tests/data b/opengate/tests/data index 1f89d1fa3f..291db9cd82 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit 1f89d1fa3f9fc593ca259ce1204c925be5a96a9f +Subproject commit 291db9cd82fedb8107ddf7b543a74d3c385a68de diff --git a/opengate/tests/src/external/pct/test101_protonct.py b/opengate/tests/src/external/pct/test101_protonct.py new file mode 100644 index 0000000000..203901a2a7 --- /dev/null +++ b/opengate/tests/src/external/pct/test101_protonct.py @@ -0,0 +1,44 @@ +import uproot +from opengate.tests import utility +from opengate.contrib.protonct.protonct import protonct + + +def compare_root_files(path1, path2, tree): + tree1 = uproot.open(path1)[tree] + df1 = tree1.arrays(library="pd") + + tree2 = uproot.open(path2)[tree] + df2 = tree2.arrays(library="pd") + + return df1.equals(df2) + + +if __name__ == "__main__": + # ===================================================== + # INITIALISATION + # ===================================================== + + paths = utility.get_default_test_paths(__file__, output_folder="test101_protonct") + + protonct(paths.output, projections=10, protons_per_projection=100, seed=1234) + + # ===================================================== + # Perform test + # ===================================================== + + path_phasespace_in = paths.output / "PhaseSpaceIn.root" + path_phasespace_out = paths.output / "PhaseSpaceOut.root" + + path_reference_phasespace_in = paths.output_ref / "PhaseSpaceIn.root" + path_reference_phasespace_out = paths.output_ref / "PhaseSpaceOut.root" + + compare_root_files(path_reference_phasespace_in, path_phasespace_in, "PhaseSpaceIn") + + is_ok = compare_root_files( + path_reference_phasespace_in, path_phasespace_in, "PhaseSpaceIn" + ) + is_ok = is_ok and compare_root_files( + path_reference_phasespace_out, path_phasespace_out, "PhaseSpaceOut" + ) + + utility.test_ok(is_ok) diff --git a/pyproject.toml b/pyproject.toml index 2b22e9b3d1..62d048853a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -60,3 +60,5 @@ phid_gammas = "opengate.bin.phid_gammas:go" phid_tac = "opengate.bin.phid_tac:go" phid_atomic_relaxation = "opengate.bin.phid_atomic_relaxation:go" phid_isomeric_transition = "opengate.bin.phid_isomeric_transition:go" + +protonct = "opengate.bin.protonct:go"