Merge pull request #608 from BishopWolf/icrp107-database

Icrp107 database

docs/requirements.txt

@@ -18,7 +18,7 @@ matplotlib
 GitPython
 colorlog
 numpy-stl
-radioactivedecay
+radioactivedecay>=0.6.1
 jsonpickle
 pandas
 requests

docs/source/user_guide/user_guide_reference_sources_generic_source.rst

@@ -343,22 +343,28 @@ normal distribution with:
 Spectra
 """""""
 
-**Discrete for gamma spectrum**
+**Discrete energy spectrum**
 
 One can configure a generic source to produce particles with energies depending on weights.
 To do so, one must provide two lists of the same size: one for energies, one for weights.
 Each energy is associated to the corresponding weight.
 Probabilities are derived from weights simply by normalizing the weights list.
 
-Several spectra are provided through the `get_rad_gamma_spectrum` function:
+1252 isotopes spectra are provided through the `get_spectrum` function:
 
 .. code:: python
 
-   spectrum = gate.sources.base.get_rad_gamma_spectrum("Lu177")
+   spectrum = gate.sources.utility.get_spectrum("Lu177", spectrum_type, database="icrp107")
 
+where ``spectrum_type`` is one of "gamma", "beta-", "beta+", "alpha", "X", "neutron",
+"auger", "IE", "alpha recoil", "anihilation", "fission", "betaD", "b-spectra". From this list,
+only b-spectra is histogram based (see next section), the rest are discrete. ``database`` can be "icrp107" or "radar".
 
-The source can be configured like this:
+ICRP107 data comes from `[ICRP, 2008. Nuclear Decay Data for Dosimetric Calculations. ICRP Publication 107. Ann. ICRP 38] <https://www.icrp.org/publication.asp?id=ICRP%20Publication%20107>`
+with the data from the `[Supplemental material] <https://journals.sagepub.com/doi/suppl/10.1177/ANIB_38_3>`.
+`[Direct link to the zipped data] <https://journals.sagepub.com/doi/suppl/10.1177/ANIB_38_3/suppl_file/P107JAICRP_38_3_Nuclear_Decay_Data_suppl_data.zip>`
 
+The source can be configured like this:
 
 .. code:: python
 
@@ -388,12 +394,12 @@ The produced particles will follow this pattern:
 One can configure a generic source to produce particles with energies according to a given histogram.
 Histograms are defined in the same way as `numpy`, using bin edges and histogram values.
 
-Several spectra are provided through the `get_rad_beta_spectrum` function.
-This data comes from `[doseinfo-radar] <https://www.doseinfo-radar.com/RADARDecay.html>`_ (`[direct link to the excel file] <https://www.doseinfo-radar.com/BetaSpec.zip>`_).
+Several spectra are provided through the `get_spectrum` function. You can use icrp107 data, or radar data.
+Radar data comes from `[doseinfo-radar] <https://www.doseinfo-radar.com/RADARDecay.html>`_ (`[direct link to the excel file] <https://www.doseinfo-radar.com/BetaSpec.zip>`_).
 
 .. code:: python
 
-   spectrum = gate.sources.base.get_rad_beta_spectrum("Lu177")
+   spectrum = gate.sources.utility.get_spectrum("Lu177", "beta-", database="radar")
 
 The source can be configured like this:
 
@@ -405,7 +411,7 @@ The source can be configured like this:
    source.energy.spectrum_energy_bin_edges = spectrum.energy_bin_edges
    source.energy.spectrum_weights = spectrum.weights
 
-For example, using this (which is what you get from `get_rad_beta_spectrum("Lu177")`):
+For example, using this (which is what you get from `get_spectrum("Lu177", "beta-", database="radar")`):
 
 .. code:: python
 

opengate/actors/base.py

@@ -305,7 +305,7 @@ def _get_docstring_for_interface(cls, output_name, interface_name):
         docstring += "\n"
         return docstring
 
-    def __init__(self, *args, **kwargs):
+    def __init__(self, *args, **kwargs) -> None:
         GateObject.__init__(self, *args, **kwargs)
         self.actor_engine = (
             None  # this is set by the actor engine during initialization

opengate/actors/digitizers.py

@@ -1083,7 +1083,7 @@ def StartSimulationAction(self):
         DigitizerBase.StartSimulationAction(self)
         DigitizerAdderActor.set_group_by_depth(self)
         if self.user_info.discretize_volume is None:
-            fatal(f'Please, set the option "discretize_volume"')
+            fatal('Please, set the option "discretize_volume"')
         depth = self.simulation.volume_manager.get_volume(
             self.discretize_volume
         ).volume_depth_in_tree

opengate/base.py

@@ -450,7 +450,7 @@ def __new__(cls, *args, **kwargs):
         new_instance = super(GateObject, cls).__new__(cls)
         return new_instance
 
-    def __init__(self, *args, simulation=None, **kwargs):
+    def __init__(self, *args, simulation=None, **kwargs) -> None:
         self._simulation = simulation
         # keep internal number of raised warnings (for debug)
         self.number_of_warnings = 0

opengate/contrib/linacs/dicomrtplan.py

@@ -17,7 +17,7 @@ def list_of_beam_sequence_ID(file):
             patient_info.value[0],
             [[("0x300C", "0x0004"), i], [("0x300A", "0x0086"), None]],
         )
-        if number_of_mu != None:
+        if number_of_mu is not None:
             id_of_beam_sequence.append(i)
     return id_of_beam_sequence
 
@@ -41,7 +41,7 @@ def read(file, cp_id="all_cp", arc_id=None):
     dose_weight_array = []
     collimation_angle_array = []
     leaf_array = []
-    if arc_id == None:
+    if arc_id is None:
         id_of_beam_sequence = list_of_beam_sequence_ID(file)
     else:
         arc_id = int(arc_id)
@@ -81,7 +81,7 @@ def read(file, cp_id="all_cp", arc_id=None):
             l_cp_id = np.arange(0, nb_cp_id, 1)
         for id in l_cp_id:
             for i, key in enumerate(key_list):
-                if data_set[i][key][id] != None:
+                if data_set[i][key][id] is not None:
                     if key == "isocenter":
                         l_parameters[i].append(np.array(data_set[i][key][id]) * mm)
                     elif (
@@ -94,8 +94,8 @@ def read(file, cp_id="all_cp", arc_id=None):
                         l_parameters[i].append(data_set[i][key][id])
                 else:
                     if id != 0:
-                        if data_set[i][key][id] == None:
-                            if l_parameters[i][id - 1] != None:
+                        if data_set[i][key][id] is None:
+                            if l_parameters[i][id - 1] is not None:
                                 l_parameters[i].append(l_parameters[i][id - 1])
                         else:
                             l_parameters[i].append(data_set[i][key][id])
@@ -120,13 +120,13 @@ def read(file, cp_id="all_cp", arc_id=None):
             leaf_blocks = [leaf_block_1, leaf_block_2]
             for i in range(len(data_set[leaves][0])):
                 for side_ID in range(2):
-                    if data_set[leaves][0][i]["leaves"][id] != None:
+                    if data_set[leaves][0][i]["leaves"][id] is not None:
                         leaf_blocks[side_ID].append(
                             data_set[leaves][side_ID][i]["leaves"][id]
                         )
                     else:
                         if id != 0:
-                            if data_set[leaves][0][i]["leaves"][id - 1] != None:
+                            if data_set[leaves][0][i]["leaves"][id - 1] is not None:
                                 leaf_blocks[side_ID].append(
                                     data_set[leaves][side_ID][i]["leaves"][id - 1]
                                 )
@@ -217,7 +217,7 @@ def read(file, cp_id="all_cp", arc_id=None):
 
 def retrieve_value_in_DICOM_RTPlan(cp, list):
     if list[0][0] in cp:
-        if list[0][1] != None:
+        if list[0][1] is not None:
             value = cp[list[0][0]].value[list[0][1]]
         else:
             value = cp[list[0][0]].value
@@ -254,7 +254,7 @@ def extract_dataset(file, beam_sequence_ID):
     data_set[dir_angle] = {"dir_angle": []}
     data_set[isocenter] = {"isocenter": []}
     data_set[dose_weights] = {"cumulative weight": []}
-    data_set[leaves] = [[{"leaves": []} for i in range(nb_leaf)] for j in range(2)]
+    data_set[leaves] = [[{"leaves": []} for _ in range(nb_leaf)] for _ in range(2)]
     data_set[MU_number] = {"MU number": 0}
     data_set[limiting_device_angle] = {"collimation angle": []}
     count = 0

opengate/contrib/phantoms/nemaiec.py

@@ -3,6 +3,7 @@
 import numpy as np
 import math
 
+from opengate import utility
 import opengate.geometry.volumes
 from opengate.utility import fatal, g4_units
 from opengate.geometry.volumes import unite_volumes
@@ -52,7 +53,7 @@ def add_iec_phantom(
     # Inside space for the water, same as the shell, with 3 mm less
     thickness = 3 * mm
     thickness_z = 10 * mm
-    interior, top_interior, c = add_iec_body(
+    interior, top_interior, _ = add_iec_body(
         simulation, f"{name}_interior", thickness, thickness_z
     )
     interior.mother = iec.name
@@ -634,7 +635,7 @@ def add_iec_phantom_vox_FIXME_TO_REMOVE(sim, name, image_filename, labels_filena
     iec = sim.add_volume("Image", name)
     iec.image = image_filename
     iec.material = "IEC_PLASTIC"
-    labels = json.loads(open(labels_filename).read())
+    labels = utility.read_json_file(labels_filename)
     iec.voxel_materials = []
     create_material(sim)
     material_list = {}
@@ -656,3 +657,32 @@ def add_iec_phantom_vox_FIXME_TO_REMOVE(sim, name, image_filename, labels_filena
         m = [labels[l]["label"], labels[l]["label"] + 1, mat]
         iec.voxel_materials.append(m)
     return iec, material_list
+
+
+def create_iec_phantom_source_vox(
+    image_filename, labels_filename, source_filename, activities=None
+):
+    if activities is None:
+        activities = {
+            "iec_sphere_10mm": 1.0,
+            "iec_sphere_13mm": 1.0,
+            "iec_sphere_17mm": 1.0,
+            "iec_sphere_22mm": 1.0,
+            "iec_sphere_28mm": 1.0,
+            "iec_sphere_37mm": 1.0,
+        }
+
+    img = itk.imread(image_filename)
+    labels = utility.read_json_file(labels_filename)
+    img_arr = itk.GetArrayViewFromImage(img)
+
+    for label in labels:
+        l = labels[label]["label"]
+        if "sphere" in label and "shell" not in label:
+            img_arr[img_arr == l] = activities[label]
+        else:
+            img_arr[img_arr == l] = 0
+
+    # The coordinate system is different from IEC analytical volume
+    # 35mm should be added in Y
+    itk.imwrite(img, source_filename)

opengate/sources/base.py

@@ -1,141 +1,8 @@
-from box import Box
-import os
-import pathlib
-import numpy as np
-import json
-
 from ..actors.base import _setter_hook_attached_to
 from ..base import GateObject, process_cls
 from ..utility import g4_units
-from ..exception import fatal
 from ..definitions import __world_name__
 
-gate_source_path = pathlib.Path(__file__).parent.resolve()
-
-# http://www.lnhb.fr/nuclear-data/module-lara/
-all_beta_plus_radionuclides = [
-    "F18",
-    "Ga68",
-    "Zr89",
-    "Na22",
-    "C11",
-    "N13",
-    "O15",
-    "Rb82",
-]
-
-
-def read_beta_plus_spectra(rad_name):
-    """
-    read the file downloaded from LNHB
-    there are 15 lines-long header to skip
-    first column is E(keV)
-    second column is dNtot/dE b+
-    WARNING : bins width is not uniform (need to scale for density)
-    """
-    filename = (
-        f"{gate_source_path}/beta_plus_spectra/{rad_name}/beta+_{rad_name}_tot.bs"
-    )
-    data = np.genfromtxt(filename, usecols=(0, 1), skip_header=15, dtype=float)
-    return data
-
-
-def compute_bins_density(bins):
-    """
-    Given a list of (energy) bins center, compute the width of each bin.
-    """
-    lower = np.roll(bins, 1)
-    lower[0] = 0
-    upper = bins
-    dx = upper - lower
-    return dx
-
-
-def get_rad_yield(rad_name):
-    if not rad_name in all_beta_plus_radionuclides:
-        return 1.0
-    data = read_beta_plus_spectra(rad_name)
-    ene = data[:, 0] / 1000  # convert from KeV to MeV
-    proba = data[:, 1]
-    cdf, total = compute_cdf_and_total_yield(proba, ene)
-    total = total * 1000  # (because was in MeV)
-    return total
-
-
-def compute_cdf_and_total_yield(data, bins):
-    """
-    Compute the CDF (Cumulative Density Function) of a list of non-uniform energy bins
-    with associated probability.
-    Also return the total probability.
-    """
-    dx = compute_bins_density(bins)
-    p = data * dx
-    total = p.sum()
-    cdf = np.cumsum(p) / total
-    return cdf, total
-
-
-def get_rad_gamma_spectrum(rad):
-    path = (
-        pathlib.Path(os.path.dirname(__file__))
-        / ".."
-        / "data"
-        / "rad_gamma_spectrum.json"
-    )
-    with open(path, "r") as f:
-        data = json.load(f)
-
-    # consider lower case
-    data = {key.lower(): value for key, value in data.items()}
-    rad = rad.lower()
-
-    if rad not in data:
-        fatal(f"get_rad_gamma_spectrum: {path} does not contain data for ion {rad}")
-
-    # select data for specific ion
-    data = Box(data[rad])
-
-    data.energies = np.array(data.energies) * g4_units.MeV
-    data.weights = np.array(data.weights)
-
-    return data
-
-
-def get_rad_beta_spectrum(rad: str):
-    path = (
-        pathlib.Path(os.path.dirname(__file__))
-        / ".."
-        / "data"
-        / "rad_beta_spectrum.json"
-    )
-    with open(path, "r") as f:
-        data = json.load(f)
-
-    if rad not in data:
-        fatal(f"get_rad_beta_spectrum: {path} does not contain data for ion {rad}")
-
-    # select data for specific ion
-    data = Box(data[rad])
-
-    bin_edges = data.energy_bin_edges
-    n = len(bin_edges) - 1
-    energies = [(bin_edges[i] + bin_edges[i + 1]) / 2 for i in range(n)]
-
-    data.energy_bin_edges = np.array(data.energy_bin_edges) * g4_units.MeV
-    data.weights = np.array(data.weights)
-    data.energies = np.array(energies)
-
-    return data
-
-
-def set_source_rad_energy_spectrum(source, rad):
-    rad_spectrum = get_rad_gamma_spectrum(rad)
-
-    source.particle = "gamma"
-    source.energy.type = "spectrum_discrete"
-    source.energy.spectrum_weights = rad_spectrum.weights
-    source.energy.spectrum_energies = rad_spectrum.energies
-
 
 class SourceBase(GateObject):
     """