Realisation of a test which read 4 randoms CP of a DICOM RT plan and …

…compare the theoretical aperture of the collimation system with the simulated one

opengate/contrib/linacs/modified_elekta_synergy_materials.db

@@ -0,0 +1,53 @@
+[Elements]
+Tungsten:   S= W   ; Z= 74. ; A= 183.84  g/mole
+Rhenium:    S= Re  ; Z= 75. ; A= 186.207 g/mole
+Copper:     S= Cu  ; Z= 29. ; A=  63.39  g/mole
+Nickel:     S= Ni  ; Z= 28. ; A=  58.693 g/mole
+Iron:       S= Fe  ; Z= 26. ; A=  55.845 g/mole
+Chromium:   S= Cr  ; Z= 24. ; A=  51.996 g/mole
+Hydrogen:   S= H   ; Z=  1. ; A=   1.01  g/mole
+Carbon:     S= C   ; Z=  6. ; A=  12.01  g/mole
+Oxygen:     S= O   ; Z=  8. ; A=  16.00  g/mole
+Aluminium:  S= Al  ; Z= 13. ; A=  26.98  g/mole
+
+[Materials]
+
+Tungsten: d=19.3 g/cm3 ; n=1 ; state=solid
+	+el: name=auto ; n=1
+
+ionizing_chamber_carbon: d=2.27 g/cm3; n=1; state=solid
+	+el: name=Carbon; n=1
+
+linac_mylar: d=1.38 g/cm3; n=3 ; state=Solid
+        +el: name=Hydrogen ; f=0.04196
+        +el: name=Carbon ; f=0.625016
+        +el: name=Oxygen; f=0.333024
+
+target_copper:	d=8.93 g/cm3; n=1; state=solid
+	+el: name=Copper ; n=1
+
+mat_leaf: d=18 g/cm3; n=3
+	+el: name=Tungsten; f=0.95
+	+el: name=Nickel;   f=0.0375
+	+el: name=Iron;     f=0.0125
+
+target_tungsten: d=19.4 g/cm3; n=2
+	+el: name=Tungsten; 	f=0.90
+	+el: name=Rhenium;	f=0.10
+
+primary_collimator_material: d=18 g/cm3; n=3
+	+el: name=Tungsten; f=0.95
+	+el: name=Nickel;   f=0.0375
+	+el: name=Iron;     f=0.0125
+
+flattening_filter_material: d=7.9 g/cm3; n=3
+	+el: name=Chromium;	f=0.17
+	+el: name=Iron;		f=0.75
+	+el: name=Nickel;	f=0.08
+
+
+linac_aluminium: d=2.7 g/cm3 ; n=1 ; state=solid
+	+el: name=Aluminium ; n=1
+
+linac_carbon: d=2.27 g/cm3 ; n=1 ; state=solid
+        +el: name=Carbon; n=1

opengate/tests/src/test069_rotation_DICOM_RT_plan.py

@@ -0,0 +1,165 @@
+import opengate as gate
+import numpy as np
+import itk
+import test069_rotation_DICOM_RT_plan_helpers as t
+from opengate.tests import utility
+
+
+def validation_test(theoretical_calculation,MC_calculation,tol =0.08):
+    print('Area from theoretical calculation for all CP (mm2): ',theoretical_calculation)
+    print('Area from MC simulations for all CP: (mm2)',MC_calculation)
+    percentage_diff = 100* (theoretical_calculation - MC_calculation)/theoretical_calculation
+    bool_percentage_diff = np.abs(percentage_diff) > tol*100
+    if np.sum(bool_percentage_diff) == 0:
+        return True
+    else :
+        return False
+
+
+def calc_MLC_aperture(x_leaf_position,y_jaws,pos_MLC=349.3,pos_jaws=470.5,SAD=1000,leaf_width = 1.85):
+
+    mm = gate.g4_units.mm
+    leaf_width = leaf_width*mm
+    left =  x_leaf_position[:80]*pos_MLC/SAD
+    right = x_leaf_position[80:]*pos_MLC/SAD
+    left[left!= 0] = left[left!= 0] - left[left!= 0]%0.5
+    right[right!=0] = right[right!=0]+ 0.5 - right[right!=0]%0.5
+
+    pos_Y_leaf = np.arange(-leaf_width*40 +leaf_width/2,leaf_width*40  -leaf_width/2 +0.01,leaf_width)
+    left[pos_Y_leaf < y_jaws[0]*pos_jaws/SAD] = 0
+    left[pos_Y_leaf > y_jaws[1] *pos_jaws/SAD] = 0
+    right[pos_Y_leaf < y_jaws[0] *pos_jaws/SAD] = 0
+    right[pos_Y_leaf > y_jaws[1]*pos_jaws/SAD] = 0
+    diff = np.array(right - left)
+
+    return np.sum(diff)*leaf_width
+
+def add_VolumeToIrradiate(sim,name,rot_volume) :
+    mm = gate.g4_units.mm
+    Box= sim.add_volume("Box", "cylinder")
+    Box.material = "G4_WATER"
+    Box.mother = name
+    Box.size = [400*mm,400*mm,400*mm]
+
+    voxel_size_x = 0.5* mm
+    voxel_size_y = 0.5 * mm
+    voxel_size_z = 400 * mm
+
+    dim_box = [400*mm/ voxel_size_x, 400*mm / voxel_size_y, 400*mm / voxel_size_z]
+    dose = sim.add_actor("DoseActor", "dose")
+    dose.output = "output/testilol.mhd"
+    dose.mother = Box.name
+    dose.size = [int(dim_box[0]),int(dim_box[1]),int(dim_box[2])]
+    dose.spacing = [voxel_size_x, voxel_size_y, voxel_size_z]
+    dose.uncertainty = False
+    dose.square = False
+    dose.hit_type = "random"
+
+    motion_tubs= sim.add_actor("MotionVolumeActor", "Move_Tubs")
+    motion_tubs.mother = Box.name
+    motion_tubs.rotations = []
+    motion_tubs.translations = []
+    for i in range(len(rot_volume)):
+        motion_tubs.rotations.append(rot_volume[i])
+        motion_tubs.translations.append([0,0,0])
+
+
+def add_alpha_source(sim,name, pos_Z,nb_part):
+    ui = sim.user_info
+    mm = gate.g4_units.mm
+    nm = gate.g4_units.nm
+    plan_source = sim.add_volume("Box","plan_alpha_source")
+    plan_source.material = "G4_Galactic"
+    plan_source.mother = name
+    plan_size = np.array([250 * mm, 148 * mm, 1 * nm])
+    plan_source.size = np.copy(plan_size)
+    plan_source.translation = [0* mm, 0 * mm, -pos_Z/2 +300*mm]
+
+
+
+    source = sim.add_source("GenericSource", "alpha_source")
+    Bq = gate.g4_units.Bq
+    MeV = gate.g4_units.MeV
+    source.particle = "alpha"
+    source.mother = plan_source.name
+    source.energy.type = "mono"
+    source.energy.mono = 1 * MeV
+    source.position.type = "box"
+    source.position.size = np.copy(plan_size)
+    source.direction.type='momentum'
+    source.direction.momentum = [0, 0, -1]
+    source.activity = nb_part * Bq / ui.number_of_threads
+
+def launch_simulation(nt,path_img,img,file,output_path,output,nb_part,src_f,vis,seg_cp,patient):
+    visu = vis
+    km = gate.g4_units.km
+    nb_cp = t.liste_CP(file)
+    nb_aleatoire = np.random.randint(0,nb_cp-1,4)
+    print('Control Points ID: ',nb_aleatoire)
+    seg_cp +=1
+    l_aperture_voxel = np.zeros(len(nb_aleatoire))
+    l_aperture_calc = np.zeros(len(nb_aleatoire))
+    for i in range(len(nb_aleatoire)):
+        cp_param = t.Dataset_DICOM_MLC_jaws(file, [nb_aleatoire[i],nb_aleatoire[i]+1], 0)
+        mean_leaves = (cp_param["Leaves"][0] + cp_param["Leaves"][1])/2
+        mean_jaws_1 = (cp_param["Y_jaws_1"][0] + cp_param["Y_jaws_1"][1])/2
+        mean_jaws_2 = (cp_param["Y_jaws_2"][0] + cp_param["Y_jaws_2"][1])/2
+        y_jaws = [mean_jaws_1,mean_jaws_2]
+        area = calc_MLC_aperture(mean_leaves,y_jaws)
+        l_aperture_calc[i] = area
+        sim = t.init_simulation(nt, cp_param,path_img, img, visu, src_f, bool_phsp=False,seg_cp=seg_cp,patient=patient)
+
+        ui = sim.user_info
+        ui.running_verbose_level = gate.logger.RUN
+
+        linac = sim.volume_manager.volumes["Linac_box"]
+        world = sim.volume_manager.volumes["world"]
+        linac.material = "G4_Galactic"
+        world.material = "G4_Galactic"
+        motion_actor = sim.get_actor_user_info("Move_LINAC")
+        rotation_volume = motion_actor.rotations
+        add_alpha_source(sim,linac.name,linac.size[2],nb_part)
+        add_VolumeToIrradiate(sim,world.name,rotation_volume)
+
+
+        dose_actor = sim.get_actor_user_info('dose')
+        dose_actor.output =  output_path / output
+
+        ui.visu = visu
+        if visu:
+            ui.visu_type = "vrml"
+        sim.physics_manager.global_production_cuts.gamma = 1 * km
+        sim.physics_manager.global_production_cuts.electron = 1 * km
+        sim.physics_manager.global_production_cuts.positron = 1 * km
+        sec = gate.g4_units.s
+        sim.run_timing_intervals = []
+        for j in range(len(rotation_volume)):
+            sim.run_timing_intervals.append([j*sec, (j+1) * sec])
+
+        sim.run(start_new_process=True)
+        path = '/home/mjacquet/Documents/Simulation_RT_plan/Simulation_MC/output/'
+        img_MC = itk.imread(output_path / output)
+        array_MC = itk.GetArrayFromImage(img_MC)
+        bool_MC = array_MC[array_MC !=0]
+        l_aperture_voxel[i] = (len(bool_MC)/4)
+    is_ok = validation_test(l_aperture_calc,l_aperture_voxel)
+    utility.test_ok(is_ok)
+
+
+if __name__ == '__main__':
+    paths = utility.get_default_test_paths(__file__)
+    output_path = paths.output
+    patient = False
+    nt = 1
+    ###### The three following variables are here to not modify the main program (the helpers) which need it ######
+    path_img = 'useless'
+    img = 'useless'
+    src_f = 'alpha'
+    ###############################################################################################################
+    output = 'img_test_069.mhd'
+    nb_part = 750000
+    seg_cp = 1
+    vis = False
+    file = '../data/DICOM_RT_plan.dcm'
+    launch_simulation(nt,path_img,img,file,output_path,output,nb_part,src_f,vis,seg_cp,patient)
+