gaga + garf standalone version numpy & torch

.github/workflows/main.yml

@@ -393,7 +393,7 @@ jobs:
             pip install torch
           fi
           pip install SimpleITK
-          pip install gaga_phsp==0.7.1
+          pip install gaga_phsp==0.7.2
           pip install dist/opengate_core-*-${PYTHONFOLDER}-${OSNAME}*_${PLATFORM}64.whl
           pip install dist/opengate-*.whl
           export GIT_SSL_NO_VERIFY=1

.gitignore

@@ -39,13 +39,11 @@ __pycache__
 /opengate/contrib/run.*
 
 /opengate/data/ensdf_230901
-
-/opengate_temporaire
+/opengate/data/isomeric_transition/save/
+/opengate/data/isomeric_transition/save2/
 
 /.pycrunch-config.yaml
-
 /opengate.egg-info/
-
 opengate_temporaire
 
 /opengate/tests/src/preload.sh
@@ -85,7 +83,6 @@ opengate_temporaire
 /opengate/tests/src/*.wrl
 /opengate/tests/src/*.gdml
 /opengate/tests/src/pb212
-/opengate/tests/src/*.gdml
 /opengate/tests/src/stl_tests
 /opengate/tests/src/label_to*
 /opengate/tests/src/test*speed*txt

docs/source/user_guide_1_intro.md

@@ -87,6 +87,6 @@ There is some additional commands lines tools that can also be used, see the [ad
 
 *Warning* they are only updated infrequently, you may have to adapt them to changes in the gate version.
 
--  [exercices](https://gitlab.in2p3.fr/davidsarrut/gate_exercices_2) (initially developed for DQPRM, French medical physics diploma)
+- [exercices](https://gitlab.in2p3.fr/davidsarrut/gate_exercices_2) (initially developed for DQPRM, French medical physics diploma)
 
 - [exercices](https://drive.google.com/drive/folders/1bcIS5OPLOBzhLo0NvrLJL5IxVQidNYCF) (initially developed for Opengate teaching)

opengate/actors/arfactors.py

@@ -42,6 +42,20 @@ def import_garf():
     return garf
 
 
+def check_channel_overlap(ch1, ch2):
+    ch1 = Box(ch1)
+    ch2 = Box(ch2)
+    if ch2.min < ch1.min < ch2.max:
+        return True
+    if ch2.min < ch1.max < ch2.max:
+        return True
+    if ch1.min < ch2.min < ch1.max:
+        return True
+    if ch1.min < ch2.max < ch1.max:
+        return True
+    return False
+
+
 class ARFTrainingDatasetActor(g4.GateARFTrainingDatasetActor, ActorBase):
     """
     The ARFTrainingDatasetActor build a root file with energy, angles, positions and energy windows
@@ -76,6 +90,17 @@ def initialize(self, simulation_engine_wr=None):
                 f"In the actor '{self.user_info.name}', the parameter 'energy_windows_actor' is {ewa.type_name}"
                 f" while it must be a DigitizerEnergyWindowsActor"
             )
+        # check overlap in channels
+        channels = ewa.channels
+        for i, ch1 in enumerate(channels):
+            for j, ch2 in enumerate(channels[i + 1 :]):
+                is_overlap = check_channel_overlap(ch1, ch2)
+                if is_overlap:
+                    fatal(
+                        f'In the actor "{self.user_info.name}", the energy channels '
+                        f"{i} and {i + 1 + j} overlap. This is not possible for ARF. \n"
+                        f"{ch1} {ch2}"
+                    )
 
     def __str__(self):
         s = f"ARFTrainingDatasetActor {self.user_info.name}"
@@ -120,20 +145,24 @@ def __init__(self, user_info):
         if self.garf is None:
             print("Cannot run GANSource")
             sys.exit()
-        # create the default detector
-        # self.user_info.arf_detector = gate.ARFDetector(self.user_info)
+
         # prepare output
         self.user_info.output_image = None
         self.g4_actor = None
         self.pth_filename = user_info.pth_filename
-        self.param = Box()
         self.nn = None
         self.model = None
         self.model_data = None
         self.batch_nb = 0
         self.detected_particles = 0
         # need a lock when the ARF is applied
         self.lock = threading.Lock()
+        # local variables
+        self.image_plane_spacing = None
+        self.image_plane_size_pixel = None
+        self.image_plane_size_mm = None
+        self.output_image = None
+        self.nb_ene = None
 
     def __str__(self):
         u = self.user_info
@@ -152,45 +181,50 @@ def __getstate__(self):
 
     def initialize(self, volume_engine=None):
         super().initialize(volume_engine)
-        # self.user_info.arf_detector.initialize(self)
         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) (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
+        self.nn, self.model = self.garf.load_nn(
+            self.pth_filename, verbose=False, gpu_mode=self.user_info.gpu_mode
+        )
+
+        # size and spacing (2D)
+        self.image_plane_spacing = np.array(
+            [self.user_info.image_spacing[0], self.user_info.image_spacing[1]]
+        )
+        self.image_plane_size_pixel = np.array(
+            [self.user_info.image_size[0], self.user_info.image_size[1]]
+        )
+        self.image_plane_size_mm = (
+            self.image_plane_size_pixel * self.image_plane_spacing
+        )
+
+        # shortcut to model_data
         self.model_data = self.nn["model_data"]
 
         # output image: nb of energy windows times nb of runs (for rotation)
-        p.nb_ene = self.model_data["n_ene_win"]
-        p.nb_runs = len(self.simulation.run_timing_intervals)
+        self.nb_ene = self.model_data["n_ene_win"]
+        nb_runs = len(self.simulation.run_timing_intervals)
         # size and spacing in 3D
-        p.image_size = [p.nb_ene, p.image_size[0], p.image_size[1]]
-        p.image_spacing = [p.image_spacing[0], p.image_spacing[1], 1]
-        # create output image as np array
-        p.output_size = [p.nb_ene * p.nb_runs, p.image_size[1], p.image_size[2]]
-        self.output_image = np.zeros(p.output_size, dtype=np.float64)
-        # compute offset
-        p.psize = [
-            p.image_size[1] * p.image_spacing[0],
-            p.image_size[2] * p.image_spacing[1],
+        self.output_image = np.array(
+            [
+                self.nb_ene,
+                self.image_plane_size_pixel[0],
+                self.image_plane_size_pixel[1],
+            ]
+        )
+        output_size = [
+            self.nb_ene * nb_runs,
+            self.output_image[1],
+            self.output_image[2],
         ]
-        p.hsize = np.divide(p.psize, 2.0)
-        p.offset = [p.image_spacing[0] / 2.0, p.image_spacing[1] / 2.0]
+        self.output_image = np.zeros(output_size, dtype=np.float64)
 
         # which device for GARF : cpu cuda mps ?
-        # we recommend CPU only
         if self.user_info.gpu_mode not in ("cpu", "gpu", "auto"):
             fatal(
                 f"the gpu_mode must be 'cpu' or 'auto' or 'gpu', while is is '{self.user_info.gpu_mode}'"
@@ -205,63 +239,57 @@ def initialize(self, volume_engine=None):
     def apply(self, actor):
         # we need a lock when the ARF is applied
         with self.lock:
-            self.apply_with_lock(actor)
+            self.arf_build_image_from_projected_points(actor)
+
+    def arf_build_image_from_projected_points(self, actor):
 
-    def apply_with_lock(self, actor):
         # get values from cpp side
         energy = np.array(actor.GetEnergy())
-        px = np.array(actor.GetPositionX())
-        py = np.array(actor.GetPositionY())
-        dx = np.array(actor.GetDirectionX())
-        dy = np.array(actor.GetDirectionY())
+        pos_x = np.array(actor.GetPositionX())
+        pos_y = np.array(actor.GetPositionY())
+        dir_x = np.array(actor.GetDirectionX())
+        dir_y = np.array(actor.GetDirectionY())
 
         # do nothing if no hits
         if energy.size == 0:
             return
 
         # convert direction in angles
-        # FIXME would it be faster on CPP side ?
         degree = g4_units.degree
-        theta = np.arccos(dy) / degree
-        phi = np.arccos(dx) / degree
+        theta = np.arccos(dir_y) / degree
+        phi = np.arccos(dir_x) / degree
 
         # update
         self.batch_nb += 1
         self.detected_particles += energy.shape[0]
 
         # build the data
-        x = np.column_stack((px, py, theta, phi, energy))
-        self.debug_nb_hits_before += len(x)
+        px = np.column_stack((pos_x, pos_y, theta, phi, energy))
+        self.debug_nb_hits_before += len(px)
 
-        # apply the neural network
+        # verbose current batch
         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 # FIXME index ?
-        w = self.garf.nn_predict(self.model, self.nn["model_data"], ax)
-
-        # positions
-        p = self.param
-        angles = x[:, 2:4]
-        t = self.garf.compute_angle_offset(angles, p.distance_to_crystal)
-        cx = x[:, 0:2]
-        cx = cx + t
-        coord = (cx + p.hsize - p.offset) / p.image_spacing[0:2]
-        coord = np.around(coord).astype(int)
-        v = coord[:, 0]
-        u = coord[:, 1]
-        u, v, w_pred = self.garf.remove_out_of_image_boundaries2(
-            u, v, w, self.user_info.image_size
+        # from projected points to image counts
+        u, v, w_pred = self.garf.arf_from_points_to_image_counts(
+            px,
+            self.model,
+            self.model_data,
+            self.user_info.distance_to_crystal,
+            self.image_plane_size_mm,
+            self.image_plane_size_pixel,
+            self.image_plane_spacing,
         )
 
         # do nothing if there is no hit in the image
         if u.shape[0] != 0:
             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)
+            s = self.nb_ene * run_id
+            img = self.output_image[s : s + self.nb_ene]
+            self.garf.image_from_coordinates_add_numpy(img, u, v, w_pred)
             self.debug_nb_hits += u.shape[0]
 
     def EndSimulationAction(self):
@@ -270,19 +298,22 @@ def EndSimulationAction(self):
         self.apply(self)
 
         # Should we keep the first slice (with all hits) ?
+        nb_slice = self.nb_ene
         if not self.user_info.enable_hit_slice:
             self.output_image = self.output_image[1:, :, :]
-            self.param.image_size[0] = self.param.image_size[0] - 1
+            # self.param.image_size[0] = self.param.image_size[0] - 1
+            nb_slice = nb_slice - 1
 
         # convert to itk image
         self.output_image = itk.image_from_array(self.output_image)
 
         # set spacing and origin like DigitizerProjectionActor
         spacing = self.user_info.image_spacing
         spacing = np.array([spacing[0], spacing[1], 1])
-        size = np.array(self.param.image_size)
-        size[0] = self.param.image_size[2]
-        size[2] = self.param.image_size[0]
+        size = np.array([0, 0, 0])
+        size[0] = self.image_plane_size_pixel[0]
+        size[1] = self.image_plane_size_pixel[1]
+        size[2] = nb_slice
         origin = -size / 2.0 * spacing + spacing / 2.0
         origin[2] = 0
         self.output_image.SetSpacing(spacing)

opengate/actors/digitizers.py

@@ -589,7 +589,10 @@ def StartSimulationAction(self):
             fatal(f"Error, the size must be 2D while it is {self.user_info.size}")
         if len(self.user_info.spacing) != 2:
             fatal(f"Error, the spacing must be 2D while it is {self.user_info.spacing}")
+        # make a copy before setting to 3 dim
+        self.user_info.size = self.user_info.size.copy()
         self.user_info.size.append(1)
+        self.user_info.spacing = self.user_info.spacing.copy()
         self.user_info.spacing.append(1)
 
         # for the moment, we cannot use this actor with several volumes

opengate/bin/opengate_tests.py

@@ -63,7 +63,8 @@ def go(test_id, random_tests):
         "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_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_3_standalone_numpy.py",  # ignored because too long (GPU)
+        "test066_spect_gaga_garf_3_standalone_torch.py",  # ignored because too long (GPU)
         "test066_spect_gaga_garf_4_analyse1.py",
         "test066_spect_gaga_garf_5_analyse2.py",
     ]

opengate/contrib/spect/siemens_intevo.py

@@ -761,7 +761,6 @@ def add_digitizer_lu177(sim, crystal_name, name):
     p1 = 112.9498 * keV
     p2 = 208.3662 * keV
     channels = [
-        {"name": "spectrum", "min": 35 * keV, "max": 588 * keV},
         *energy_windows_peak_scatter("peak113", "scatter1", "scatter2", p1, 0.2, 0.1),
         *energy_windows_peak_scatter("peak208", "scatter3", "scatter4", p2, 0.2, 0.1),
     ]
@@ -820,7 +819,6 @@ def add_digitizer_tc99m(sim, crystal_name, name):
     # energy windows (Energy range. 35-588 keV)
     cc = digitizer.add_module("DigitizerEnergyWindowsActor", f"{name}_energy_window")
     channels = [
-        {"name": f"spectrum", "min": 3 * keV, "max": 160 * keV},
         {"name": f"scatter", "min": 108.57749938965 * keV, "max": 129.5924987793 * keV},
         {"name": f"peak140", "min": 129.5924987793 * keV, "max": 150.60751342773 * keV},
     ]

opengate/engines.py

@@ -2,15 +2,11 @@
 import random
 import sys
 import os
-
-# from multiprocessing import Process, set_start_method, Manager
-# import queue
 import weakref
 from box import Box
 from anytree import PreOrderIter
 
 import opengate_core as g4
-
 from .exception import fatal, warning
 from .decorators import requires_fatal, requires_warning
 from .logger import log
@@ -525,11 +521,10 @@ def get_actor(self, name):
         return self.actors[name]
 
     def create_actors(self):
-        for (
-            ui
-        ) in (
-            self.simulation_engine_wr().simulation.actor_manager.user_info_actors.values()
-        ):
+        # consider the priority value of the actors
+        uia = self.simulation_engine_wr().simulation.actor_manager.user_info_actors
+        sorted_uia = sorted(uia.values(), key=lambda d: d.priority)
+        for ui in sorted_uia:
             actor = new_element(ui, self.simulation_engine_wr().simulation)
             log.debug(f"Actor: initialize [{ui.type_name}] {ui.name}")
             actor.initialize(self.simulation_engine_wr)