From b75e96595115305546a1c1471222f1e0c095e5d7 Mon Sep 17 00:00:00 2001 From: Martina Favaretto Date: Mon, 23 Jan 2023 16:59:00 +0100 Subject: [PATCH] created tests for tps source Treatment plan source file was missing from the branch. Added. small modifications. Still trying to figure out the error in spot positions Fixed bug with TPS rotation and created tests Removed old version of the helpers. Now the function that reads the treatment plan txt file returns also the gantry angle Created tests for: optics in air, absolute dose in water and gantry rotation. Working version of TPSource Removed full MA polinomia from the beamline model added test for weights and comparison pbs-tps exposed G4UserLimits and StepLimiter Corrected bugs in TPSource and created working tests for optics, absolute dose, ranges and weights Improved test metrics and TPS readability test optics in vbl Refactored TPS init and renamed beamline properties changed TPS initialization minor changes, mainly estetics' turn off visualization reduced the size of the detectors Further reduced the detector's side Included energy spread in the TPS model and examples Corrected inconsistency in test for absolute dose. Added function to test dose grid shape and spacing to helpers_test accounted for correction in dose grid (PR#137) accounted for PR#137 corrected scale_dose function for case image size 1,1,x multiprocessing: do we need a timeout in p.join()? spot weighted on the beamset number of particles, not on the beam ones. created separate neme dictionary for particles to use in user_limit Changed queue initialization for multiprocessing added check to avoid error when beam spot has almost zero weight removed step limiter from branch modified threshold --- opengate/SimulationEngine.py | 37 +- opengate/__init__.py | 3 + opengate/helpers_beamline.py | 69 ++ opengate/helpers_rt_plan.py | 698 ++++++++++++++++++ opengate/helpers_tests.py | 255 ++++++- opengate/source/TreatmentPlanSource.py | 135 ++++ .../tests/src/test059_TPSource_abs_dose.py | 194 +++++ .../tests/src/test059_TPSource_gantry_rot.py | 201 +++++ opengate/tests/src/test059_TPSource_optics.py | 225 ++++++ .../tests/src/test059_TPSource_optics_vbl.py | 211 ++++++ .../tests/src/test059_TPSource_range_ref.py | 161 ++++ .../tests/src/test059_TPSource_weights.py | 206 ++++++ 12 files changed, 2386 insertions(+), 9 deletions(-) create mode 100644 opengate/helpers_beamline.py create mode 100644 opengate/helpers_rt_plan.py create mode 100644 opengate/source/TreatmentPlanSource.py create mode 100644 opengate/tests/src/test059_TPSource_abs_dose.py create mode 100644 opengate/tests/src/test059_TPSource_gantry_rot.py create mode 100644 opengate/tests/src/test059_TPSource_optics.py create mode 100644 opengate/tests/src/test059_TPSource_optics_vbl.py create mode 100644 opengate/tests/src/test059_TPSource_range_ref.py create mode 100644 opengate/tests/src/test059_TPSource_weights.py diff --git a/opengate/SimulationEngine.py b/opengate/SimulationEngine.py index 1170a0ecf..7e3878e7b 100644 --- a/opengate/SimulationEngine.py +++ b/opengate/SimulationEngine.py @@ -4,7 +4,14 @@ import sys import os from .ExceptionHandler import * -from multiprocessing import Process, set_start_method, Queue +from multiprocessing import ( + Process, + set_start_method, + Manager, + Queue, + active_children, + cpu_count, +) from opengate_core import G4RunManagerFactory from .Decorators import requires_fatal from .helpers import fatal, warning @@ -131,12 +138,28 @@ def start(self): # (the "force" option is needed for notebooks) set_start_method("fork", force=True) # set_start_method("spawn") - q = Queue() + q = Manager().Queue() + # q = Queue() p = Process(target=self.init_and_start, args=(q,)) + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {q.full()}") + print("---start process---") p.start() + import time + + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {q.full()}") + while len(active_children()) >= cpu_count() + 4: + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + time.sleep(0.01) + print(q.full()) self.state = "started" - p.join() + p.join() # (timeout=10) # timeout might be needed self.state = "after" + print("AFTER") output = q.get() else: output = self.init_and_start(None) @@ -184,7 +207,15 @@ def init_and_start(self, queue): output.store_hook_log(self) output.current_random_seed = self.current_random_seed if queue is not None: + print("--- in process, before put ---") + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {queue.full()}") queue.put(output) + print("--- in process, after put ---") + print(f"Active children: {len(active_children())}") + print(f"CPU count: {cpu_count()}") + print(f"Queue full: {queue.full()}") return None else: return output diff --git a/opengate/__init__.py b/opengate/__init__.py index 6bdb1e077..78efc3abc 100644 --- a/opengate/__init__.py +++ b/opengate/__init__.py @@ -10,6 +10,8 @@ from .helpers_tests import * from .helpers_tests_root import * from .helpers_transform import * +from .helpers_beamline import * +from .helpers_rt_plan import * # main mechanism for the 'elements': source, actor, volume from .UserInfo import * @@ -17,6 +19,7 @@ from .source.SourceBase import * from .actor.ActorBase import * from .geometry.VolumeBase import * +from .source.TreatmentPlanSource import * # main object from .Simulation import * diff --git a/opengate/helpers_beamline.py b/opengate/helpers_beamline.py new file mode 100644 index 000000000..02995878a --- /dev/null +++ b/opengate/helpers_beamline.py @@ -0,0 +1,69 @@ +# N.B: distances in mm, degrees in rad + + +class Rashi: + pass + + +class RangeMod: + pass + + +class BeamlineModel: + def __init__(self): + self.name = None + self.radiation_types = [] + self.rm = None # range modulator + self.rashi = None + # Nozzle entrance to Isocenter distance + self.distance_nozzle_iso = 0 # mm + # SMX (X bending magnet) to Isocenter distance + self.distance_stearmag_to_isocenter_x = float( + "inf" + ) # default infinity for parallel beams + # SMY (Y bending magnet) to Isocenter distance + self.distance_stearmag_to_isocenter_y = float("inf") + # polinomial coefficients + self.energy_mean_coeffs = [0] + self.energy_spread_coeffs = [0] + self.sigma_x_coeffs = [0] + self.theta_x_coeffs = [0] + self.epsilon_x_coeffs = [0] + self.sigma_y_coeffs = [0] + self.theta_y_coeffs = [0] + self.epsilon_y_coeffs = [0] + # convergence + self.conv_x = 0 + self.conv_y = 0 + + def _polynomial_map(self, base, coeff): + # coeff are given with decreasing degree (coeff[0]->max degree) + polyDegree = len(coeff) + exp = list(range(polyDegree)) + exp.reverse() + + return sum([c * (base ** (i)) for c, i in zip(coeff, exp)]) + + def get_energy(self, nominal_energy): + return self._polynomial_map(nominal_energy, self.energy_mean_coeffs) + + def get_sigma_energy(self, nominal_energy): + return self._polynomial_map(nominal_energy, self.energy_spread_coeffs) + + def get_sigma_x(self, energy): + return self._polynomial_map(energy, self.sigma_x_coeffs) + + def get_theta_x(self, energy): + return self._polynomial_map(energy, self.theta_x_coeffs) + + def get_epsilon_x(self, energy): + return self._polynomial_map(energy, self.epsilon_x_coeffs) + + def get_sigma_y(self, energy): + return self._polynomial_map(energy, self.sigma_y_coeffs) + + def get_theta_y(self, energy): + return self._polynomial_map(energy, self.theta_y_coeffs) + + def get_epsilon_y(self, energy): + return self._polynomial_map(energy, self.epsilon_y_coeffs) diff --git a/opengate/helpers_rt_plan.py b/opengate/helpers_rt_plan.py new file mode 100644 index 000000000..0d8ba9492 --- /dev/null +++ b/opengate/helpers_rt_plan.py @@ -0,0 +1,698 @@ +#!/usr/bin/env python3 +# ----------------------------------------------------------------------------- +# Copyright (C): MedAustron GmbH, ACMIT Gmbh and Medical University Vienna +# This software is distributed under the terms +# of the GNU Lesser General Public Licence (LGPL) +# See LICENSE for further details +# ----------------------------------------------------------------------------- + +import pydicom +import os +import re +import numpy as np + +# from utils.dose_info import dose_info +import logging + +logger = logging.getLogger(__name__) + + +def is_close(x, y, eps=1e-6): + sumabs = np.abs(x) + np.abs(y) + absdif = np.abs(x - y) + ok = (sumabs == 0) or (absdif < eps * 0.5 * sumabs) + return ok + + +class spot_info(object): + def __init__(self, xiec, yiec, w, e): + self.xiec = xiec + self.yiec = yiec + self.w = w + self.energy = e + self.particle_name = None + self.beamFraction = None + self.t0 = None + self.t1 = None + + def get_msw(self, t0, t1): + return self.w + + +class layer_info(object): + def __init__(self, ctrlpnt, j, cumsumchk=[], verbose=False, keep0=False): + self._cp = ctrlpnt + if verbose: + logger.debug("{}. control point with type {}".format(j, type(self._cp))) + for k in self._cp.keys(): + if pydicom.datadict.dictionary_has_tag(k): + kw = pydicom.datadict.keyword_for_tag(k) + else: + kw = "(UNKNOWN)" + logger.debug("k={} keyword={}".format(k, kw)) + nspot = int(self._cp.NumberOfScanSpotPositions) + # assert(self._cp.NominalBeamEnergyUnit == 'MEV') + if nspot == 1: + self.w = np.array([float(self._cp.ScanSpotMetersetWeights)]) + else: + self.w = np.array([float(w) for w in self._cp.ScanSpotMetersetWeights]) + assert nspot == len(self.w) + assert nspot * 2 == len(self._cp.ScanSpotPositionMap) + # self.cpindex = int(self._cp.ControlPointIndex) + # self.spotID = str(self._cp.ScanSpotTuneID) + cmsw = float(self._cp.CumulativeMetersetWeight) + if cumsumchk: + logger.debug( + "CumulativeMetersetWeight={0:14.8g} sum of previous spots={1:14.8g} diff={2:14.8g}".format( + cmsw, cumsumchk[0], cmsw - cumsumchk[0] + ) + ) + assert is_close(cmsw, cumsumchk[0]) + # self.npainting = int(self._cp.NumberOfPaintings) + xy = np.array([float(pos) for pos in self._cp.ScanSpotPositionMap]).reshape( + nspot, 2 + ) + self.x = np.array(xy[:, 0]) + self.y = np.array(xy[:, 1]) + if not keep0: + mask = self.w > 0.0 + self.w = self.w[mask] + self.x = self.x[mask] + self.y = self.y[mask] + + wsum = np.sum(self.w) + logger.debug( + "layer number {} has {} spots, of which {} with zero weight, cumsum={}, sum(w)={}".format( + j, len(self.w), np.sum(self.w <= 0), cmsw, wsum + ) + ) + + cumsumchk[0] += wsum + + @property + def energy(self): + # DICOM specifies energy per nucleon, Gate wants total kinetic energy + return float(self._cp.NominalBeamEnergy) + + @property + def tuneID(self): + return str(self._cp.ScanSpotTuneID) + + @property + def npainting(self): + return int(self._cp.NumberOfPaintings) + + @property + def mswtot(self): + return np.sum(self.w) + + @property + def nspots(self): + return len(self.w) + + @property + def weights(self): + return self.w + + @property + def spots(self): + e = self.energy + return [spot_info(x, y, w, e) for (x, y, w) in zip(self.x, self.y, self.w)] + + def get_spots(self, t0=None, t1=None): + e = self.energy + return [spot_info(x, y, w, e) for (x, y, w) in zip(self.x, self.y, self.w)] + + +class beam_info(object): + # def __init__(self,beam,rd,i,keep0=False): + def __init__(self, beam, i, override_number, keep0=False): + logger.debug("loading {}th beam".format(i)) + self._dcmbeam = beam # the DICOM beam object + self._warnings = list() # will hopefully stay empty + logger.debug("trying to access first control point") + self._icp0 = beam.IonControlPointSequence[0] # convenience: first control point + self._beam_number_is_fishy = ( + override_number # workaround for buggy TPSs, e.g. PDM + ) + self._index = i # the index in the beam sequence + self._layers = list() + mswchk = self.FinalCumulativeMetersetWeight + cumsumchk = [0.0] + logger.debug("going to read all layers") + for j, icp in enumerate(self._dcmbeam.IonControlPointSequence): + li = layer_info(icp, j, cumsumchk, False, keep0) + if 0.0 < li.mswtot or keep0: + self._layers.append(li) + logger.debug("survived reading all layers") + if not is_close(mswchk, cumsumchk[0]): + raise ValueError( + "final cumulative msw {} != sum of spot msw {}".format( + mswchk, cumsumchk[0] + ) + ) + logger.debug("survived cumulative MSW check") + + def GetAndClearWarnings(self): + # return and clear + w = self._warnings[:] + self._warnings = list() + return w + + @property + def FinalCumulativeMetersetWeight(self): + return float(self._dcmbeam.FinalCumulativeMetersetWeight) + + @property + def PatientSupportAngle(self): + return float(self._icp0.PatientSupportAngle) + + @property + def patient_angle(self): + return float(self._icp0.PatientSupportAngle) + + @property + def IsoCenter(self): + if "IsocenterPosition" in self._icp0: + if len(self._icp0.IsocenterPosition) == 3.0: + return [float(xyz) for xyz in self._icp0.IsocenterPosition] + else: + msg = "Got corrupted isocenter = '{}'; assuming [0,0,0] for now, keep fingers crossed.".format( + self._icp0.IsocenterPosition + ) + else: + msg = "No isocenter specified in treatment plan; assuming [0,0,0] for now, keep fingers crossed." + logger.error(msg) + self._warnings.append(msg) + # FIXME: what to do else? Cry? Throw segfaults? Drink bad coffee? + return [0.0, 0.0, 0.0] + + @property + def Name(self): + # TODO: the Name and name properties are identical, keep only one of them. + return str(self._dcmbeam.BeamName) + + @property + def Number(self): + # TODO: the Number and number properties are identical, keep only one of them. + nr = ( + str(self._index + 1) + if self._beam_number_is_fishy + else str(self._dcmbeam.BeamNumber) + ) + return nr + + @property + def name(self): + # TODO: the Name and name properties are identical, keep only one of them. + return str(self._dcmbeam.BeamName) + + @property + def number(self): + # TODO: the Number and number properties are identical, keep only one of them. + nr = ( + str(self._index + 1) + if self._beam_number_is_fishy + else str(self._dcmbeam.BeamNumber) + ) + return nr + + @property + def RadiationType(self): + radtype = str(self._dcmbeam.RadiationType) + radtypeOpengate = None + if radtype == "ION": + ionZ = str(self._dcmbeam.RadiationAtomicNumber) + ionA = str(self._dcmbeam.RadiationMassNumber) + ionQ = str(self._dcmbeam.RadiationChargeState) + radtype = "_".join(["ION", ionZ, ionA, ionQ]) + radtypeOpengate = f"ion {ionZ} {ionA}" + return [radtype, radtypeOpengate] + + @property + def gantry_angle(self): + return float(self._icp0.GantryAngle) + + @property + def TreatmentMachineName(self): + if "TreatmentMachineName" in self._dcmbeam: + return str(self._dcmbeam.TreatmentMachineName) + # RayStation 8b exports anonymized treatment plans without the treatment machine name! + if np.isclose(self.gantry_angle, 0.0): + # FIXME: should be solved in a way that works for any clinic, not just MedAustron + ducktape = str("IR2VBL") + elif np.isclose(self.gantry_angle, 90.0): + # FIXME: should be solved in a way that works for any clinic, not just MedAustron + ducktape = str("IR2HBL") + else: + raise ValueError( + "treatment machine name is missing and gantry angle {} does not enable a good guess".format( + self.gantry_angle + ) + ) + msg = "treatment machine name for beam name={} number={} missing in treatment plan, guessing '{}' from gantry angle {}".format( + self.name, self.number, ducktape, self.gantry_angle + ) + logger.error(msg) + self._warnings.append(msg) + return ducktape # ugly workaround! FIXME! + + @property + def SnoutID(self): + if "SnoutSequence" in self._dcmbeam: + return str(self._dcmbeam.SnoutSequence[0].SnoutID) + # FIXME: what to do else? + return str("NA") + + @property + def SnoutPosition(self): + if "SnoutPosition" in self._dcmbeam: + return float(self._icp0.SnoutPosition) + # FIXME: what to do else? + return str("NA") + + @property + def NumberOfRangeModulators(self): + return int(self._dcmbeam.NumberOfRangeModulators) + + @property + def RangeModulatorIDs(self): + if self.NumberOfRangeModulators > 0: + return [rm.RangeModulatorID for rm in self._dcmbeam.RangeModulatorSequence] + return list() + + @property + def NumberOfRangeShifters(self): + return int(self._dcmbeam.NumberOfRangeShifters) + + @property + def RangeShifterIDs(self): + if self.NumberOfRangeShifters > 0: + return [str(rs.RangeShifterID) for rs in self._dcmbeam.RangeShifterSequence] + return list() + + @property + def NumberOfEnergies(self): + return len( + set( + [icp.NominalBeamEnergy for icp in self._dcmbeam.IonControlPointSequence] + ) + ) + + @property + def nlayers(self): + return len(self._layers) + + @property + def layers(self): + return self._layers + + @property + def nspots(self): + return sum([l.nspots for l in self.layers]) + + @property + def mswtot(self): + return sum([l.mswtot for l in self._layers]) + + @property + def PrimaryDosimeterUnit(self): + return str(self._dcmbeam.PrimaryDosimeterUnit) + + +def sequence_check(obj, attr, nmin=1, nmax=0, name="object"): + logger.debug("checking that {} has attribute {}".format(name, attr)) + assert hasattr(obj, attr) + seq = getattr(obj, attr) + logger.debug( + "{} has length {}, will check if it >={} and <={}".format( + name, len(seq), nmin, nmax + ) + ) + assert len(seq) >= nmin + assert nmax == 0 or len(seq) <= nmax + + +class beamset_info(object): + """ + This class reads a DICOM 'RT Ion Plan Storage' file and collects related information such as TPS dose files. + It does NOT (yet) try to read a reffered structure set and/or CT images. + This acts as a wrapper (all DICOM access on the plan file happens here). This has a few advantages over direct + DICOM access in the other modules: + * we can deal with different "DICOM dialects" here; some TPSs may store their plans in different ways. + * if 'private tags' need to be taken into account then we can also do that here. + * We can make a similar class, with the same attributes, for a treatment plan stored in a different format, e.g. for research, commissioning or QA purposes. + + Then the rest of the code can work with that in the same way. + """ + + patient_attrs = ["Patient ID", "Patient Name", "Patient Birth Date", "Patient Sex"] + plan_req_attrs = [ + "RT Plan Label", + "SOP Instance UID", + "Referring Physician Name", + "Plan Intent", + ] + plan_opt_attrs = ["Operators Name", "Reviewer Name", "Review Date", "Review Time"] + plan_attrs = plan_req_attrs + plan_opt_attrs + bs_attrs = [ + "Number Of Beams", + "RT Plan Label", + "Prescription Dose", + "Target ROI Name", + "Radiation Type", + "Treatment Machine(s)", + ] + + def __init__(self, rpfp): + self._warnings = list() # will hopefully stay empty + self._beam_numbers_corrupt = False # e.g. PDM does not define beam numbers + self._rp = pydicom.read_file(rpfp) + self._rpfp = rpfp + logger.debug("beamset: survived reading DICOM file {}".format(rpfp)) + self._rpdir = os.path.dirname(rpfp) + self._rpuid = str(self._rp.SOPInstanceUID) + self._dose_roiname = ( + None # stays None for CT-less plans, e.g. commissioning plans + ) + self._dose_roinumber = ( + None # stays None for CT-less plans, e.g. commissioning plans + ) + logger.debug("beamset: going to do some checks") + self._chkrp() + logger.debug("beamset: survived check, loading beams") + self._beams = [ + beam_info(b, i, self._beam_numbers_corrupt) + for i, b in enumerate(self._rp.IonBeamSequence) + ] + logger.debug("beamset: DONE") + + def GetAndClearWarnings(self): + # return a copy + for b in self._beams: + # bwarnings = b.GetAndClearWarnings() + for w in b.GetAndClearWarnings(): + if w not in self._warnings: + self._warnings.append(w) + allw = self._warnings[:] + self._warnings = list() + return allw + + def __getitem__(self, k): + if type(k) == int: + if k >= 0 and k < len(self._beams): + return self._beams[k] + else: + raise IndexError( + "attempt to get nonexisting beam with index {}".format(k) + ) + for b in self._beams: + if str(k) == b.name or str(k) == b.number: + return b + raise KeyError("attempt to get nonexisting beam with key {}".format(k)) + + def _chkrp(self): + if "SOPClassUID" not in self._rp: + raise IOError("bad DICOM file {},\nmissing SOPClassUID".format(self._rpfp)) + sop_class_name = pydicom.uid.UID_dictionary[self._rp.SOPClassUID][0] + if sop_class_name != "RT Ion Plan Storage": + raise IOError( + "bad plan file {},\nwrong SOPClassUID: {}='{}',\nexpecting an 'RT Ion Plan Storage' file instead.".format( + self._rpfp, self._rp.SOPClassUID, sop_class_name + ) + ) + missing_attrs = list() + for a in ["IonBeamSequence"] + self.plan_req_attrs + self.patient_attrs: + b = a.replace(" ", "") + if not hasattr(self._rp, b): + missing_attrs.append(b) + if missing_attrs: + raise IOError( + "bad plan file {},\nmissing keys: {}".format( + self._rpfp, ", ".join(missing_attrs) + ) + ) + # self._get_rds() + if hasattr(self._rp, "DoseReferenceSequence"): + sequence_check(self._rp, "DoseReferenceSequence", 1, 1) + if hasattr(self._rp.DoseReferenceSequence[0], "ReferencedROINumber"): + self._dose_roinumber = int( + self._rp.DoseReferenceSequence[0].ReferencedROINumber + ) + if self._dose_roinumber is None: + logger.info( + "no target ROI specified (probably because of missing DoseReferenceSequence)" + ) + sequence_check(self._rp, "IonBeamSequence", 1, 0) + sequence_check(self._rp, "FractionGroupSequence", 1, 1) + frac0 = self._rp.FractionGroupSequence[0] + sequence_check( + frac0, + "ReferencedBeamSequence", + len(self._rp.IonBeamSequence), + len(self._rp.IonBeamSequence), + ) + number_set = set() + for dcmbeam in self._rp.IonBeamSequence: + nr = int(dcmbeam.BeamNumber) + if nr < 0: + self._beam_numbers_corrupt = True + logger.error( + "CORRUPT INPUT: found a negative beam number {}.".format(nr) + ) + if nr in number_set: + self._beam_numbers_corrupt = True + logger.error( + "CORRUPT INPUT: found same beam number {} for multiple beams.".format( + nr + ) + ) + number_set.add(nr) + if self._beam_numbers_corrupt: + msg = "Beam numbers are corrupt! Will override them with simple enumeration, keep fingers crossed." + logger.error(msg) + self._warnings.append(msg) + logger.debug("checked planfile, looks like all attributes are available...") + + @property + def mswtot(self): + return sum([b.mswtot for b in self._beams]) + + @property + def name(self): + # It looks like RTPlanLabel is for the beamset, + # and RTPlanName is for the entire plan including possibly several beamsets. + # the RTPlanName is not exported anymore in RayStation 8b, so let's forget about the plan name + # some anonymization methods are butchering all useful labeling information, even the label and name of an RT plan. + return str(self._rp.get("RTPlanLabel", "anonymized")) + + @property + def fields(self): + # GATE synomym for 'beams' + return self._beams + + @property + def beams(self): + return self._beams + + @property + def uid(self): + return self._rpuid + + @property + def beam_angles(self): + return [str(bm.gantry_angle) for bm in self._beams] + + @property + def beam_names(self): + return [str(bm.Name) for bm in self._beams] + + @property + def beam_numbers(self): + return [str(bm.Number) for bm in self._beams] + + @property + def Nfractions(self): + # FIXME: some evil DICOM plan files have "NumberOfFractionsPlanned" equal to zero. Force this to be one, or is this zero somehow meaningful & useful? + nfrac = int(self._rp.FractionGroupSequence[0].NumberOfFractionsPlanned) + if nfrac > 0: + return nfrac + logger.error("Got Nfractions={} ???! Using nfrac=1 instead.".format(nfrac)) + return 1 + + @property + def target_ROI_name(self): + return self._dose_roiname + + @target_ROI_name.setter + def target_ROI_name(self, roiname): + self._dose_roiname = roiname + + @property + def target_ROI_number(self): + return self._dose_roinumber + + @property + def prescription_dose(self): + if hasattr(self._rp, "DoseReferenceSequence"): + if hasattr(self._rp.DoseReferenceSequence[0], "TargetPrescriptionDose"): + return float(self._rp.DoseReferenceSequence[0].TargetPrescriptionDose) + return "NA" + + @property + def plan_label(self): + return str(self._rp.get("RTPlanLabel", "anonymized")) + + @property + def sanitized_plan_label(self): + badchars = re.compile("[^a-zA-Z0-9_]") + return re.sub(badchars, "_", self.plan_label) + + @property + def patient_info(self): + return dict( + [ + (a, str(getattr(self._rp, a.replace(" ", "")))) + for a in self.patient_attrs + ] + ) + + @property + def plan_info(self): + reqs = dict( + [ + (a, str(getattr(self._rp, a.replace(" ", "")))) + for a in self.plan_req_attrs + ] + ) + opts = dict( + [ + ( + a, + "NA" + if not hasattr(self._rp, a.replace(" ", "")) + else str(getattr(self._rp, a.replace(" ", ""))), + ) + for a in self.plan_opt_attrs + ] + ) + reqs.update(opts) + return reqs + + @property + def bs_info(self): + info = dict( + [ + ("Number Of Beams", str(len(self._beams))), + ("RT Plan Label", self.plan_label), + ("Prescription Dose", str(self.prescription_dose)), + ("Target ROI Name", str(self.target_ROI_name)), + ( + "Radiation Type", + ", ".join( + set([str(beam.RadiationType[0]) for beam in self._beams]) + ), + ), + ( + "Radiation Type Opengate", + ", ".join( + set([str(beam.RadiationType[1]) for beam in self._beams]) + ), + ), + ( + "Treatment Machine(s)", + ", ".join( + set([str(beam.TreatmentMachineName) for beam in self._beams]) + ), + ), + ] + ) + dirty = self.plan_label + sanitized = self.sanitized_plan_label + if dirty != sanitized: + info["Sanitized RT Plan Label"] = sanitized + return info + + def __repr__(self): + s = "\nPLAN\n\t" + "\n\t".join( + ["{0:30s}: {1}".format(a, self.plan_info[a]) for a in self.plan_attrs] + ) + s += "\nPATIENT\n\t" + "\n\t".join( + ["{0:30s}: {1}".format(a, self.patient_info[a]) for a in self.patient_attrs] + ) + s += "\nBEAMSET\n\t" + "\n\t".join( + ["{0:30s}: {1}".format(a, self.bs_info[a]) for a in self.bs_attrs] + ) + return s + + +def spots_info_from_txt(txtFile, ionType): + # initialize empty variables + nFields = 0 + ntot = [] + energies = [] + nSpots = [] + spots = [] + start_index = [] + G = 0 + + # read content + with open(txtFile, "r") as f: + lines = f.readlines() + + # get plan's info + for i, line in enumerate(lines): + if line.startswith("###GantryAngle"): + l = lines[i + 1].split("\n")[0] + G = int(l) + if line.startswith("##NumberOfFields"): + l = lines[i + 1].split("\n")[0] + nFields = int(l) + if line.startswith("###FinalCumulativeMeterSetWeight"): + l = lines[i + 1].split("\n")[0] + ntot.append(float(l)) + if line.startswith("####Energy"): + l = lines[i + 1].split("\n")[0] + energies.append(float(l)) + if line.startswith("####NbOfScannedSpots"): + l = lines[i + 1].split("\n")[0] + nSpots.append(int(l)) + if line.startswith("####X Y Weight"): + start_index.append(i + 1) + + np = sum(ntot) + for k in range(nFields): + # np = ntot[k] + for i in range(len(energies)): + e = energies[i] + print(f"ENERGY: {e}") + start = start_index[i] + end = start_index[i] + nSpots[i] + for j in range(start, end): + l = lines[j].split("\n")[0].split() + spot = spot_info(float(l[0]), float(l[1]), float(l[2]), e) + spot.beamFraction = float(l[2]) / np + spot.particle_name = ionType + spots.append(spot) + + return spots, np, energies, G + + +def get_spots_from_beamset(beamset): + rad_type = beamset.bs_info["Radiation Type Opengate"] + spots_array = [] + mswtot = beamset.mswtot + for beam in beamset.beams: + # mswtot = beam.mswtot + for energy_layer in beam.layers: + for spot in energy_layer.spots: + nPlannedSpot = spot.w + spot.beamFraction = ( + nPlannedSpot / mswtot + ) # nr particles planned for the spot/tot particles planned for the beam + spot.particle_name = rad_type + spots_array.append(spot) + return spots_array + + +# vim: set et softtabstop=4 sw=4 smartindent: diff --git a/opengate/helpers_tests.py b/opengate/helpers_tests.py index ed9bc540d..46ef0ceb9 100644 --- a/opengate/helpers_tests.py +++ b/opengate/helpers_tests.py @@ -209,7 +209,7 @@ def plot_img_y(ax, img, label): data = itk.GetArrayViewFromImage(img) y = np.sum(data, 2) y = np.sum(y, 0) - x = np.arange(len(y)) * img.GetSpacing()[2] + x = np.arange(len(y)) * img.GetSpacing()[1] ax.plot(x, y, label=label) ax.legend() @@ -219,7 +219,7 @@ def plot_img_x(ax, img, label): data = itk.GetArrayViewFromImage(img) y = np.sum(data, 1) y = np.sum(y, 0) - x = np.arange(len(y)) * img.GetSpacing()[2] + x = np.arange(len(y)) * img.GetSpacing()[0] ax.plot(x, y, label=label) ax.legend() @@ -975,10 +975,13 @@ def plot_gauss_fit(positionVec, dose, fit, show=False): return fig -def create_position_vector(length, spacing): +def create_position_vector(length, spacing, centered=True): # cretae position vector, with origin in the image plane's center width = length * spacing - positionVec = np.arange(0, width, spacing) - width / 2 + spacing / 2 + if centered: + positionVec = np.arange(0, width, spacing) - width / 2 + spacing / 2 + else: + positionVec = np.arange(0, width, spacing) return positionVec @@ -1007,14 +1010,31 @@ def read_mhd(filename): return data, spacing, shape -def create_2D_Edep_colorMap(filepath, show=False): +def plot2D(twodarray, label, show=False): + fig = plt.figure(figsize=(20, 20)) + ax = fig.add_subplot(111) + ax.set_title(label) + plt.imshow(twodarray) + ax.set_aspect("equal") + plt.colorbar(orientation="vertical") + if show: + plt.show() + return fig + + +def create_2D_Edep_colorMap(filepath, show=False, axis="z"): img = itk.imread(str(filepath)) data = itk.GetArrayViewFromImage(img) fig = plt.figure(figsize=(20, 20)) ax = fig.add_subplot(111) ax.set_title("colorMap") - plt.imshow(data[0, :, :]) + if axis == "z": + plt.imshow(data[0, :, :]) + elif axis == "x": + plt.imshow(data[:, :, 0]) + else: + plt.imshow(data[:, 0, :]) ax.set_aspect("equal") plt.colorbar(orientation="vertical") if show: @@ -1192,6 +1212,229 @@ def test_weights(expected_ratio, mhd_1, mhd_2, thresh=0.1): return is_ok +def test_tps_spot_size_positions(data, ref, spacing, thresh=0.1, abs_tol=0.3): + if not np.array_equal(data.size, ref.size): + print("Images do not have the same size") + return False + ok = True + # beam along x + size = data.shape + # get gaussian fit of Edep only around the i-th spot + param_y_out, _, _ = gate.extract_gauss_param_1D(data, size[1], spacing[1], axis=0) + param_y_ref, _, _ = gate.extract_gauss_param_1D(ref, size[1], spacing[1], axis=0) + + param_z_out, _, _ = gate.extract_gauss_param_1D(data, size[0], spacing[2], axis=1) + param_z_ref, _, _ = gate.extract_gauss_param_1D(ref, size[0], spacing[2], axis=1) + + # check positions + print("Check position of the spot") + print(f" opengate: ({param_y_out[1]:.2f},{param_z_out[1]:.2f})") + print(f" gate: ({param_y_ref[1]:.2f},{param_z_ref[1]:.2f})") + + diffmY = param_y_out[1] - param_y_ref[1] # / param_y_ref[1] + diffmZ = param_z_out[1] - param_z_ref[1] # / param_z_ref[1] + mean_diff = np.mean([diffmY, diffmZ]) + + if ( + (abs(diffmY) > 2 * abs_tol) + or (abs(diffmZ) > 2 * abs_tol) + or (abs(mean_diff) > abs_tol) + ): + print( + f"\033[91m Position error above threshold. DiffX={diffmY:.2f}, diffY={diffmZ:.2f}, threshold is 0.3mm \033[0m" + ) + ok = False + + # check sizes + print("Check size of the spot") + print(f" opengate: ({param_y_out[2]:.2f},{param_z_out[2]:.2f})") + print(f" gate: ({param_y_ref[2]:.2f},{param_z_ref[2]:.2f})") + + diffsY = (param_y_out[2] - param_y_ref[2]) / param_y_ref[2] + diffsZ = (param_z_out[2] - param_z_ref[2]) / param_z_ref[2] + + if (diffsY > thresh) or (diffsZ > thresh): + print("\033[91m Size error above threshold \033[0m") + ok = False + + return ok + + +def scale_dose(path, scaling, outpath): + img_mhd_in = itk.imread(path) + data = itk.GetArrayViewFromImage(img_mhd_in) + dose = data * scaling + spacing = img_mhd_in.GetSpacing() + img = gate.itk_image_view_from_array(dose) + img.SetSpacing(spacing) + itk.imwrite(img, outpath) + return outpath + + +def check_dose_grid_geometry(dose_mhd_path, dose_actor): + img = itk.imread(dose_mhd_path) + data = itk.GetArrayViewFromImage(img) + shape = data.shape + spacing = img.GetSpacing() + shape_ref = tuple(np.flip(dose_actor.size)) + spacing_ref = dose_actor.spacing + + ok = True + if shape != shape_ref: + print(f"{shape=} not the same as {shape_ref=}!") + ok = False + + if spacing != spacing_ref: + print(f"{spacing=} not the same as {spacing_ref=}!") + ok = False + + return ok + + +def arangeDx(dx, xV, includeUB=False, lb=[], ub=[]): + if not lb: + lb = np.amin(xV) + if not ub: + ub = np.amax(xV) + if includeUB: + x_int = np.arange(lb, ub + dx / 10, dx) + else: + x_int = np.arange(lb, ub, dx) + return x_int + + +def interpolate1Dprofile(xV, dV, dx=0.01, interpolMethod="cubic"): + f = scipy.interpolate.interp1d( + xV, dV, kind=interpolMethod, fill_value="extrapolate" + ) + xVfine = arangeDx(dx, xV, includeUB=True, lb=np.amin(xV), ub=np.amax(xV)) + dVfine = f(xVfine) + return xVfine, dVfine + + +def getRange(xV, dV, percentLevel=0.8): + dx = 0.01 + xVfine, dVfine = interpolate1Dprofile(xV, dV, dx, "cubic") + + indMaxFine = np.argmax(dVfine) + indR80 = np.argmax( + np.logical_and( + xVfine > xVfine[indMaxFine], dVfine <= percentLevel * dVfine[indMaxFine] + ) + ) + r80 = xVfine[indR80] + dAtR80 = dVfine[indR80] + + return (r80, dAtR80) + + +def get_range_from_image(volume, shape, spacing, axis="y"): + x1, d1 = get_1D_profile(volume, shape, spacing, axis=axis) + r, _ = getRange(x1, d1) + + return r + + +def compareRange( + volume1, + volume2, + shape1, + shape2, + spacing1, + spacing2, + axis1="y", + axis2="y", + thresh=2.0, +): + ok = True + x1, d1 = get_1D_profile(volume1, shape1, spacing1, axis=axis1) + x2, d2 = get_1D_profile(volume2, shape2, spacing2, axis=axis2) + + print("---RANGE80---") + r1, _ = getRange(x1, d1) + r2, _ = getRange(x2, d2) + print(r1) + print(r2) + diff = abs(r2 - r1) + + if diff > thresh: + print(f"\033[91mRange difference is {diff}mm, threshold is {thresh}mm \033[0m") + ok = False + + return ok + + +def get_1D_profile(data, shape, spacing, axis="z"): + if axis == "x": + d1 = np.sum(np.sum(data, 1), 0) + x1 = create_position_vector(shape[2], spacing[0], centered=False) + + if axis == "y": + d1 = np.sum(np.sum(data, 2), 0) + x1 = create_position_vector(shape[1], spacing[1], centered=False) + + if axis == "z": + d1 = np.sum(np.sum(data, 2), 1) + x1 = create_position_vector(shape[0], spacing[2], centered=False) + + return x1, d1 + + +def compare_dose_at_points( + pointsV, + dose1, + dose2, + shape1, + shape2, + spacing1, + spacing2, + axis1="z", + axis2="z", + rel_tol=0.03, +): + ok = True + s1 = 0 + s2 = 0 + x1, doseV1 = get_1D_profile(dose1, shape1, spacing1, axis=axis1) + x2, doseV2 = get_1D_profile(dose2, shape2, spacing2, axis=axis2) + # plt.plot(x1, doseV1) + # plt.plot(x2, doseV2) + # plt.show() + for p in pointsV: + # get dose at the position p [mm] + cp1 = min(x1, key=lambda x: abs(x - p)) + d1_p = doseV1[np.where(x1 == cp1)] + + cp2 = min(x2, key=lambda x: abs(x - p)) + d2_p = doseV2[np.where(x2 == cp2)] + + s1 += d1_p + s2 += d2_p + + print(abs(s1 - s2) / s2) + + # print(f"Dose difference at {p} mm is {diff_pc}%") + if abs(s1 - s2) / s2 > rel_tol: + print(f"\033[91mDose difference above threshold \033[0m") + ok = False + return ok + + +def assert_img_sum(img1, img2, sum_tolerance=5): + data1 = itk.GetArrayViewFromImage(img1).ravel() + data2 = itk.GetArrayViewFromImage(img2).ravel() + + s1 = np.sum(data1) + s2 = np.sum(data2) + if s1 == 0 and s2 == 0: + t = 0 + else: + t = np.fabs((s1 - s2) / s1) * 100 + b = t < sum_tolerance + print_test(b, f"Img sums {s1} vs {s2} : {t:.2f} % (tol {sum_tolerance:.2f} %)") + return b + + def check_diff(value1, value2, tolerance, txt): diff = np.fabs(value1 - value2) / value1 * 100 t = diff < tolerance diff --git a/opengate/source/TreatmentPlanSource.py b/opengate/source/TreatmentPlanSource.py new file mode 100644 index 000000000..6089fcc58 --- /dev/null +++ b/opengate/source/TreatmentPlanSource.py @@ -0,0 +1,135 @@ +import numpy as np +from scipy.spatial.transform import Rotation +import opengate as gate + + +class TreatmentPlanSource: + def __init__(self, name, sim): + self.name = name + # self.mother = None + self.rotation = Rotation.identity() + self.translation = [0, 0, 0] + self.spots = None + self.beamline_model = None + self.n_sim = 0 + self.sim = sim # simulation obj to which we want to add the tp source + + def __del__(self): + pass + + def set_particles_to_simulate(self, n_sim): + self.n_sim = n_sim + + def set_spots(self, spots): + self.spots = spots + + def set_spots_from_rtplan(self, rt_plan_path): + beamset = gate.beamset_info(rt_plan_path) + gantry_angle = beamset.beam_angles[0] + spots = gate.get_spots_from_beamset(beamset) + self.spots = spots + self.rotation = Rotation.from_euler("z", gantry_angle, degrees=True) + + def set_beamline_model(self, beamline): + self.beamline_model = beamline + + def initialize_tpsource(self): + # some alias + spots_array = self.spots + sim = self.sim + nSim = self.n_sim + beamline = self.beamline_model + self.d_nozzle_to_iso = beamline.distance_nozzle_iso + self.d_stearMag_to_iso_x = beamline.distance_stearmag_to_isocenter_x + self.d_stearMag_to_iso_y = beamline.distance_stearmag_to_isocenter_y + + # mapping factors between iso center plane and nozzle plane (due to steering magnets) + cal_proportion_factor = ( + lambda d_magnet_iso: 1 + if (d_magnet_iso == float("inf")) + else (d_magnet_iso - self.d_nozzle_to_iso) / d_magnet_iso + ) + self.proportion_factor_x = cal_proportion_factor(self.d_stearMag_to_iso_x) + self.proportion_factor_y = cal_proportion_factor(self.d_stearMag_to_iso_y) + + # initialize a pencil beam for each spot + for i, spot in enumerate(spots_array): + nspot = np.round( + spot.beamFraction * nSim + ) # simulate a fraction of the beam particles for this spot + if nspot == 0: + continue + + source = sim.add_source("PencilBeamSource", f"{self.name}_spot_{i}") + + # set energy + source.energy.type = "gauss" + source.energy.mono = beamline.get_energy(nominal_energy=spot.energy) + source.energy.sigma_gauss = beamline.get_sigma_energy( + nominal_energy=spot.energy + ) + + source.particle = spot.particle_name + source.position.type = "disc" # pos = Beam, shape = circle + sigma + + # # set mother + # if self.mother is not None: + # source.mother = self.mother + + # POSITION: + source.position.translation = self._get_pbs_position(spot) + + # ROTATION: + source.position.rotation = self._get_pbs_rotation(spot) + + # add weight + # source.weight = -1 + source.n = nspot + + # set optics parameters + source.direction.partPhSp_x = [ + beamline.get_sigma_x(spot.energy), + beamline.get_theta_x(spot.energy), + beamline.get_epsilon_x(spot.energy), + beamline.conv_x, + ] + source.direction.partPhSp_y = [ + beamline.get_sigma_y(spot.energy), + beamline.get_theta_y(spot.energy), + beamline.get_epsilon_y(spot.energy), + beamline.conv_y, + ] + + def _get_pbs_position(self, spot): + # (x,y) referr to isocenter plane. + # Need to be corrected to referr to nozzle plane + pos = [ + (spot.xiec) * self.proportion_factor_x, + (spot.yiec) * self.proportion_factor_y, + self.d_nozzle_to_iso, + ] + # Gantry angle = 0 -> source comes from +y and is positioned along negative side of y-axis + # https://opengate.readthedocs.io/en/latest/source_and_particle_management.html + + position = (self.rotation * Rotation.from_euler("x", np.pi / 2)).apply( + pos + ) + self.translation + + return position + + def _get_pbs_rotation(self, spot): + # by default the source points in direction z+. + # Need to account for SM direction deviation and rotation thoward isocenter (270 deg around x) + # then rotate of gantry angle + rotation = [0.0, 0.0, 0.0] + beta = np.arctan(spot.yiec / self.d_stearMag_to_iso_y) + alpha = np.arctan(spot.xiec / self.d_stearMag_to_iso_x) + rotation[0] = -np.pi / 2 + beta + rotation[2] = -alpha + + # apply gantry angle + spot_rotation = ( + self.rotation * Rotation.from_euler("xyz", rotation) + ).as_matrix() + + return spot_rotation diff --git a/opengate/tests/src/test059_TPSource_abs_dose.py b/opengate/tests/src/test059_TPSource_abs_dose.py new file mode 100644 index 000000000..89a319e23 --- /dev/null +++ b/opengate/tests/src/test059_TPSource_abs_dose.py @@ -0,0 +1,194 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.translation = [1148 * mm, 0.0, 0.0] +box.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [500 * mm, 500 * mm, 400 * mm] +phantom.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() +phantom.translation = [-200.0, 0.0, 0] +phantom.material = "G4_WATER" +phantom.color = [0, 0, 1, 1] + +# roos chamber +roos = sim.add_volume("Tubs", "roos") +roos.mother = phantom.name +roos.material = "G4_WATER" +roos.rmax = 7.8 +roos.rmin = 0 +roos.dz = 200 +roos.color = [1, 0, 1, 1] + + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" #'QGSP_BIC_HP_EMZ' #"FTFP_INCLXX_EMZ" + +sim.set_cut("world", "all", 1000 * km) + + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "abs_dose_roos.mhd" +dose.mother = roos.name +dose.size = [1, 1, 800] +dose.spacing = [15.6, 15.6, 0.5] +dose.hit_type = "random" +dose.gray = True + + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2HBL = gate.BeamlineModel() +IR2HBL.name = None +IR2HBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2HBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +IR2HBL.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +IR2HBL.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +IR2HBL.sigma_x_coeffs = [2.3335753978880014] +IR2HBL.theta_x_coeffs = [0.0002944903217664001] +IR2HBL.epsilon_x_coeffs = [0.0007872786903040108] +IR2HBL.sigma_y_coeffs = [1.9643343053823967] +IR2HBL.theta_y_coeffs = [0.0007911780133478402] +IR2HBL.epsilon_y_coeffs = [0.0024916149017600447] + +## --------START PENCIL BEAM SCANNING---------- ## +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) + +nSim = 200000 # 328935 # particles to simulate per beam + +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-F5x5cm_E120MeVn.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2HBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = gate.scale_dose( + str(dose.output).replace(".mhd", "_dose.mhd"), + ntot / nSim, + output_path / "threeDdoseWaternew.mhd", +) + +# ABSOLUTE DOSE + +# read output and ref +img_mhd_out = itk.imread(dose_path) +img_mhd_ref = itk.imread( + ref_path / "idc-PHANTOM-roos-F5x5cm_E120MeVn-PLAN-Physical.mhd" +) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +shape = data.shape +spacing = img_mhd_out.GetSpacing() +spacing_ref = np.flip(img_mhd_ref.GetSpacing()) + +ok = gate.assert_img_sum( + img_mhd_out, + img_mhd_ref, +) + +points = 400 - np.linspace(10, 14, 9) +ok = ( + gate.compare_dose_at_points( + points, + data, + data_ref, + shape, + data_ref.shape, + spacing, + spacing_ref, + axis1="z", + axis2="x", + rel_tol=0.03, + ) + and ok +) + +# 1D +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_out, "x profile", axis="z") +# gate.plot_img_axis(ax, img_mhd_ref, "x ref", axis="z") +# plt.show() + +# fig.savefig(output_path / "dose_profiles_water.png") + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_gantry_rot.py b/opengate/tests/src/test059_TPSource_gantry_rot.py new file mode 100644 index 000000000..ec3d7f595 --- /dev/null +++ b/opengate/tests/src/test059_TPSource_gantry_rot.py @@ -0,0 +1,201 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +## ---------- DEFINE BEAMLINE MODEL -------------## +beamline = gate.BeamlineModel() +beamline.name = None +beamline.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +beamline.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +beamline.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +beamline.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +beamline.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +beamline.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +beamline.sigma_x_coeffs = [2.3335753978880014] +beamline.theta_x_coeffs = [0.0002944903217664001] +beamline.epsilon_x_coeffs = [0.0007872786903040108] +beamline.sigma_y_coeffs = [1.9643343053823967] +beamline.theta_y_coeffs = [0.0007911780133478402] +beamline.epsilon_y_coeffs = [0.0024916149017600447] + +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) +nSim = 20000 # particles to simulate per beam +# rt_plan = ref_path / "RP1.2.752.243.1.1.20220406175810679.4500.52008_tagman.dcm" +# beamset = gate.beamset_info(rt_plan) +# G = float(beamset.beam_angles[0]) + + +## ---- VBL Nozzle --- +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.rotation = Rotation.from_euler("x", -90, degrees=True).as_matrix() +box.translation = [0.0, -1148 * mm, 0.0] +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# Rashi +rashi = sim.add_volume("Box", "rashi") +rashi.mother = box.name +rashi.size = [500 * mm, 500 * mm, 5 * mm] +rashi.translation = [0.0, 0.0, 200 * mm] +rashi.material = "G4_LUCITE" +rashi.color = [1, 0, 1, 1] + +## ---- HBL Nozzle --- +box_rot = sim.add_volume("Box", "box_rot") +gate.copy_user_info(box, box_rot) +box_rot.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box_rot.translation = [1148.0, 0.0, 1000.0] + +nozzle_rot = sim.add_volume("Box", "nozzle_rot") +gate.copy_user_info(nozzle, nozzle_rot) +nozzle_rot.mother = box_rot.name + +rashi_rot = sim.add_volume("Box", "rashi_rot") +gate.copy_user_info(rashi, rashi_rot) +rashi_rot.mother = box_rot.name + +# ----------------------------------- + +# target 1 VBL +phantom = sim.add_volume("Box", "phantom") +phantom.size = [324 * mm, 324 * mm, 324 * mm] +phantom.translation = [0 * mm, 0.0, 0.0] +phantom.material = "G4_WATER" +phantom.color = [0, 0, 1, 1] + +# target 2 HBL +phantom_rot = sim.add_volume("Box", "phantom_rot") +gate.copy_user_info(phantom, phantom_rot) +phantom_rot.rotation = Rotation.from_euler("z", 90, degrees=True).as_matrix() +phantom_rot.translation = [0.0, 0.0, 1000.0] + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "testTPSgantry.mhd" +dose.mother = phantom.name +dose.size = [162, 1620, 162] +dose.spacing = [2.0, 0.2, 2.0] +dose.hit_type = "random" +dose.gray = True + +dose_rot = sim.add_actor("DoseActor", "doseInXYZ_rot") +gate.copy_user_info(dose, dose_rot) +dose_rot.mother = phantom_rot.name +dose_rot.output = output_path / "testTPSganry_rot.mhd" + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" +sim.set_cut("world", "all", 1000 * km) + +# add TPSources +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-1D_HBL_120.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("VBL", sim) +tps.set_beamline_model(beamline) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", 0, degrees=True) +tps.initialize_tpsource() + +tps_rot = gate.TreatmentPlanSource("HBL", sim) +tps_rot.set_beamline_model(beamline) +tps_rot.set_particles_to_simulate(nSim) +tps_rot.set_spots(spots) +tps_rot.rotation = Rotation.from_euler("z", G, degrees=True) +tps_rot.translation = [0.0, 0.0, 1000.0] +tps_rot.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## + +# ABSOLUTE DOSE +# ok = gate.assert_images( +# dose.output, +# dose_rot.output, +# stat, +# tolerance=50, +# ignore_value=0, +# ) +ok = True + +# read output and ref +img_mhd_out = itk.imread(dose_rot.output) +img_mhd_ref = itk.imread(dose.output) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +spacing = img_mhd_out.GetSpacing() + +# Range 80 +ok = ( + gate.compareRange(data, data_ref, data.shape, data_ref.shape, spacing, spacing) + and ok +) + +# 1D plots +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_out, "y profile", axis="y") +# gate.plot_img_axis(ax, img_mhd_ref, "y ref", axis="y") +# fig.savefig(output_path / "dose_profiles_water.png") +# plt.show() + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_optics.py b/opengate/tests/src/test059_TPSource_optics.py new file mode 100644 index 000000000..be2152c28 --- /dev/null +++ b/opengate/tests/src/test059_TPSource_optics.py @@ -0,0 +1,225 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.translation = [1148 * mm, 0.0, 0.0] +box.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [300 * mm, 310 * mm, 310 * mm] +phantom.translation = [150 * mm, 0.0, 0.0] +# phantom.rotation = Rotation.from_euler('x',-90,degrees=True).as_matrix() +phantom.material = "G4_AIR" +phantom.color = [0, 0, 1, 1] + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" # "Shielding_EMZ"#"FTFP_INCLXX_EMZ" + +# p.physics_list_name = "QGSP_BIC_EMZ" +sim.set_cut("world", "all", 1000 * km) + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "testTPSoptics.mhd" +dose.mother = phantom.name +dose.size = [30, 620, 620] +dose.spacing = [10.0, 0.5, 0.5] +dose.hit_type = "random" +dose.gray = True + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2HBL = gate.BeamlineModel() +IR2HBL.name = None +IR2HBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2HBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +IR2HBL.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +IR2HBL.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +IR2HBL.sigma_x_coeffs = [-0.00011142901344618727, 2.346946879501544] +IR2HBL.theta_x_coeffs = [-3.6368814874049214e-07, 0.0003381328996152591] +IR2HBL.epsilon_x_coeffs = [3.1292233857396716e-06, 0.0004117718840152502] +IR2HBL.sigma_y_coeffs = [-0.0004009682717802152, 2.0124504979960225] +IR2HBL.theta_y_coeffs = [-8.437400716390318e-07, 0.000892426821944524] +IR2HBL.epsilon_y_coeffs = [-8.757558864087579e-08, 0.00250212397239695] + +## --------START PENCIL BEAM SCANNING---------- ## +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) +# nSim = 328935 # particles to simulate per beam +nSim = 20000 +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-gate_test51_TP_1_old.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2HBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +# rt_plan = ref_path / "RP1.2.752.243.1.1.20230119115736709.2000.75541.dcm" +# tps.set_spots_from_rtplan(rt_plan) # no need to set rotation here +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = gate.scale_dose( + str(dose.output).replace(".mhd", "_dose.mhd"), + ntot / nSim, + output_path / "threeDdoseAirSpots.mhd", +) + +# SPOT POSITIONS COMPARISON +# read output and ref +img_mhd_out = itk.imread(dose_path) +img_mhd_ref = itk.imread( + ref_path / "idc-PHANTOM-air_box-gate_test51_TP_1-PLAN-Physical.mhd" +) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +shape = data.shape +spacing = img_mhd_out.GetSpacing() + +# spot comparison (N.B x and z are inverted in np array!) +# spots in the plan file +yz = [ + 50, + -50, + -50, + -50, + -100, + 0, + 100, + 0, + 50, + 50, + 0, + 100, + -100, + -100, + 0, + -100, + 100, + -100, + 0, + 0, + -50, + 50, + -100, + 100, + 100, + 100, +] + +yzM = np.array(yz).reshape(int(len(yz) / 2), 2) +# convert from mm (wrt image center) to voxel +spot_y = [int(y / dose.spacing[1]) + int(dose.size[1] / 2) for y in yzM[:, 0]] +spot_z = [int(z / dose.spacing[1]) + int(dose.size[1] / 2) for z in yzM[:, 1]] + +thresh = 0.8 ## OSS: need step limiter + +# # 1D +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_out, "z profile", axis="z") +# #gate.plot_img_axis(ax, img_mhd_out, "x profile", axis="x") +# gate.plot_img_axis(ax, img_mhd_out, "y profile", axis="y") + + +# # fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +# gate.plot_img_axis(ax, img_mhd_ref, "z ref", axis="z") +# #gate.plot_img_axis(ax, img_mhd_ref, "x ref", axis="x") +# gate.plot_img_axis(ax, img_mhd_ref, "y ref", axis="y") + +# fig.savefig(output_path / "dose_profiles_spots.png") + + +ok = True +for i in range(1, shape[2], shape[2] // 3): + # check i-th slab + print(f"Airslab nr. {i}") + # gate.plot2D(data[:, :, i], "2D Edep opengate", show=True) + # gate.plot2D(data_ref[:, :, i], "2D Edep gate", show=True) + for y, z in zip(spot_y, spot_z): + # i = 0 + print(f" ({y:.2f},{z:.2f})") + # 'cut' the slab around the spot expected in y,z + w = 30 # cut window's half size + d_out = data[z - w : z + w, y - w : y + w, i : i + 1] + d_ref = data_ref[z - w : z + w, y - w : y + w, i : i + 1] + ok = ( + gate.test_tps_spot_size_positions( + d_out, d_ref, spacing, thresh=thresh, abs_tol=1.0 + ) + and ok + ) + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_optics_vbl.py b/opengate/tests/src/test059_TPSource_optics_vbl.py new file mode 100644 index 000000000..3a58ea53e --- /dev/null +++ b/opengate/tests/src/test059_TPSource_optics_vbl.py @@ -0,0 +1,211 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.rotation = Rotation.from_euler("x", -90, degrees=True).as_matrix() +box.translation = [0.0, -1148 * mm, 0.0] +box.material = "Air" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [300 * mm, 310 * mm, 310 * mm] +phantom.rotation = Rotation.from_euler("z", -90, degrees=True).as_matrix() +phantom.translation = [0.0, -150 * mm, 0.0] +phantom.material = "G4_AIR" +phantom.color = [0, 0, 1, 1] + + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" +# p.physics_list_name = "QGSP_BIC_EMZ" +sim.set_cut("world", "all", 1000 * km) + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "TPS_optics_vbl.mhd" +dose.mother = phantom.name +dose.size = [30, 620, 620] +dose.spacing = [10.0, 0.5, 0.5] +dose.hit_type = "random" +dose.gray = True + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2VBL = gate.BeamlineModel() +IR2VBL.name = None +IR2VBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2VBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2VBL.distance_stearmag_to_isocenter_x = 6700.00 +# SMY to Isocenter distance +IR2VBL.distance_stearmag_to_isocenter_y = 7420.00 +# polinomial coefficients +IR2VBL.energy_mean_coeffs = [11.91893485094217, -9.539517997860457] +IR2VBL.energy_spread_coeffs = [0.0004790681841295621, 5.253257865904452] +IR2VBL.sigma_x_coeffs = [-0.00011142901344618727, 2.346946879501544] +IR2VBL.theta_x_coeffs = [-3.6368814874049214e-07, 0.0003381328996152591] +IR2VBL.epsilon_x_coeffs = [3.1292233857396716e-06, 0.0004117718840152502] +IR2VBL.sigma_y_coeffs = [-0.0004009682717802152, 2.0124504979960225] +IR2VBL.theta_y_coeffs = [-8.437400716390318e-07, 0.000892426821944524] +IR2VBL.epsilon_y_coeffs = [-8.757558864087579e-08, 0.00250212397239695] + + +## --------START PENCIL BEAM SCANNING---------- ## +# nSim = 328935 # particles to simulate per beam +nSim = 20000 +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan4Gate-gate_test51tps_v.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2VBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = gate.scale_dose( + str(dose.output).replace(".mhd", "_dose.mhd"), + ntot / nSim, + output_path / "threeDdoseAirSpots_vbl.mhd", +) + +# SPOT POSITIONS COMPARISON +# read output and ref +img_mhd_out = itk.imread(dose_path) +img_mhd_ref = itk.imread( + ref_path / "idc-PHANTOM-air_box_vbl-gate_test51tps_v-PLAN-Physical.mhd" +) +data = itk.GetArrayViewFromImage(img_mhd_out) +data_ref = itk.GetArrayViewFromImage(img_mhd_ref) +shape = data.shape +spacing = img_mhd_out.GetSpacing() + + +# spot comparison (N.B x and z are inverted in np array!) +# spots in the plan file +yz = [ + 0, + 50, + 0, + 0, + -62.507, + -21.669, + -61.951, + -72.23, + 50, + -50, + 50, + 0, + 50, + 50, + -50, + 50, +] + +yzM = np.array(yz).reshape(int(len(yz) / 2), 2) +# convert from mm (wrt image center) to voxel +spot_y = [int(y / dose.spacing[1]) + int(dose.size[1] / 2) for y in yzM[:, 0]] +spot_z = [int(z / dose.spacing[1]) + int(dose.size[1] / 2) for z in yzM[:, 1]] + +thresh = 0.8 ## OSS: need step limiter + +# 1D +fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +gate.plot_img_axis(ax, img_mhd_out, "z profile", axis="z") +gate.plot_img_axis(ax, img_mhd_out, "x profile", axis="x") +gate.plot_img_axis(ax, img_mhd_out, "y profile", axis="y") + +# fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) +gate.plot_img_axis(ax, img_mhd_ref, "z ref", axis="z") +gate.plot_img_axis(ax, img_mhd_ref, "x ref", axis="x") +gate.plot_img_axis(ax, img_mhd_ref, "y ref", axis="y") +fig.savefig(output_path / "dose_profiles_spots_vbl.png") + +ok = True +for i in range(1, shape[2], shape[2] // 3): + # check i-th slab + print(f"Airslab nr. {i}") + # gate.plot2D(data[:, :, i], "2D Edep opengate", show=True) + # gate.plot2D(data_ref[:, :, i], "2D Edep gate", show=True) + for y, z in zip(spot_y, spot_z): + # i = 0 + print(f" ({y:.2f},{z:.2f})") + # 'cut' the slab around the spot expected in y,z + w = 30 # cut window's half size + d_out = data[z - w : z + w, y - w : y + w, i : i + 1] + d_ref = data_ref[z - w : z + w, y - w : y + w, i : i + 1] + ok = ( + gate.test_tps_spot_size_positions( + d_out, d_ref, spacing, thresh=thresh, abs_tol=1.0 + ) + and ok + ) + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_range_ref.py b/opengate/tests/src/test059_TPSource_range_ref.py new file mode 100644 index 000000000..ef4d641eb --- /dev/null +++ b/opengate/tests/src/test059_TPSource_range_ref.py @@ -0,0 +1,161 @@ +import opengate as gate +import itk +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +## ------ INITIALIZE SIMULATION ENVIRONMENT ---------- ## +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 12365478910 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +rad = gate.g4_units("rad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +# nozzle box +box = sim.add_volume("Box", "box") +box.size = [500 * mm, 500 * mm, 1000 * mm] +box.translation = [1148 * mm, 0.0, 0.0] +box.rotation = Rotation.from_euler("y", -90, degrees=True).as_matrix() +box.material = "Vacuum" +box.color = [0, 0, 1, 1] + +# nozzle WET +nozzle = sim.add_volume("Box", "nozzle") +nozzle.mother = box.name +nozzle.size = [500 * mm, 500 * mm, 2 * mm] +nozzle.material = "G4_WATER" + +# target +phantom = sim.add_volume("Box", "phantom") +phantom.size = [500 * mm, 500 * mm, 400 * mm] +phantom.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() +phantom.translation = [-200.0, 0.0, 0] +phantom.material = "G4_WATER" +phantom.color = [0, 0, 1, 1] + +# roos chamber +peak_finder = sim.add_volume("Tubs", "peak_finder") +peak_finder.mother = phantom.name +peak_finder.material = "G4_WATER" +peak_finder.rmax = 40.3 +peak_finder.rmin = 0 +peak_finder.dz = 200 +peak_finder.color = [1, 0, 1, 1] + + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" #'QGSP_BIC_HP_EMZ' #"FTFP_INCLXX_EMZ" +sim.set_cut("world", "all", 1000 * km) + + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInXYZ") +dose.output = output_path / "dose_peak_finder.mhd" +dose.mother = peak_finder.name +dose.size = [1, 1, 8000] +dose.spacing = [80.6, 80.6, 0.05] +dose.hit_type = "random" +dose.gray = True + + +## ---------- DEFINE BEAMLINE MODEL -------------## +IR2HBL = gate.BeamlineModel() +IR2HBL.name = None +IR2HBL.radiation_types = "ion 6 12" +# Nozzle entrance to Isocenter distance +IR2HBL.distance_nozzle_iso = 1300.00 # 1648 * mm#1300 * mm +# SMX to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_x = 2000.00 +# SMY to Isocenter distance +IR2HBL.distance_stearmag_to_isocenter_y = 2000.00 +# polinomial coefficients +IR2HBL.energy_mean_coeffs = [12.0, -9.54] +IR2HBL.energy_spread_coeffs = [0.00048, 5.2532] +IR2HBL.sigma_x_coeffs = [2.33] +IR2HBL.theta_x_coeffs = [0.00029] +IR2HBL.epsilon_x_coeffs = [0.00078] +IR2HBL.sigma_y_coeffs = [1.96] +IR2HBL.theta_y_coeffs = [0.00079] +IR2HBL.epsilon_y_coeffs = [0.0024] + +## --------START PENCIL BEAM SCANNING---------- ## +# NOTE: HBL means that the beam is coming from -x (90 degree rot around y) +nSim = 20000 # 328935 # particles to simulate per beam +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "PlanCentralSpot_1440MeV.txt", "ion 6 12" +) +tps = gate.TreatmentPlanSource("RT_plan", sim) +tps.set_beamline_model(IR2HBL) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("z", G, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True +# start simulation +sim.run() +output = sim.output + +## -------------END SCANNING------------- ## +# print results at the end +stat = output.get_actor("Stats") +print(stat) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +## ------ TESTS -------## +dose_path = str(dose.output).replace(".mhd", "_dose.mhd") + +# RANGE + +# read output and ref +img_mhd_out = itk.imread(dose_path) +data = itk.GetArrayViewFromImage(img_mhd_out) +shape = data.shape +spacing = img_mhd_out.GetSpacing() + +# Range 80 +range80_gate9_E120MeV = 367.06 +range_opengate = gate.get_range_from_image(data, data.shape, spacing, axis="z") + +thresh = 2.0 * mm +ok = True +if abs(range_opengate - range80_gate9_E120MeV) > thresh: + ok = False + + +gate.test_ok(ok) diff --git a/opengate/tests/src/test059_TPSource_weights.py b/opengate/tests/src/test059_TPSource_weights.py new file mode 100644 index 000000000..6f655abcf --- /dev/null +++ b/opengate/tests/src/test059_TPSource_weights.py @@ -0,0 +1,206 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from scipy.spatial.transform import Rotation +import os + +paths = gate.get_default_test_paths(__file__, "gate_test044_pbs") +output_path = paths.output / "output_test051_rtp" +ref_path = paths.output_ref / "test051_ref" + +# create the simulation +sim = gate.Simulation() + +# main options +ui = sim.user_info +ui.g4_verbose = False +ui.g4_verbose_level = 1 +ui.visu = False +ui.random_seed = 123654789 +ui.random_engine = "MersenneTwister" + +# units +km = gate.g4_units("km") +cm = gate.g4_units("cm") +mm = gate.g4_units("mm") +um = gate.g4_units("um") +MeV = gate.g4_units("MeV") +Bq = gate.g4_units("Bq") +nm = gate.g4_units("nm") +deg = gate.g4_units("deg") +mrad = gate.g4_units("mrad") + +# add a material database +sim.add_material_database(paths.gate_data / "HFMaterials2014.db") + +# change world size +world = sim.world +world.size = [600 * cm, 500 * cm, 500 * cm] + +## FIRST DETECTOR ## +# box +# translation and rotation like in the Gate macro +box1 = sim.add_volume("Box", "box1") +box1.size = [10 * cm, 10 * cm, 105 * cm] +box1.translation = [0 * cm, 0 * cm, 52.5 * cm] +box1.material = "Vacuum" +box1.color = [0, 0, 1, 1] + +# phantoms +m = Rotation.identity().as_matrix() + +phantom = sim.add_volume("Box", "phantom_a_1") +phantom.mother = "box1" +phantom.size = [100 * mm, 100 * mm, 50 * mm] +phantom.translation = [0 * mm, 0 * mm, -500 * mm] +phantom.rotation = m +phantom.material = "G4_AIR" +phantom.color = [1, 0, 1, 1] + + +# add dose actor +dose = sim.add_actor("DoseActor", "doseInYZ_1") +filename = "phantom_a_1.mhd" +dose.output = output_path / filename +dose.mother = "phantom_a_1" +dose.size = [250, 250, 1] +dose.spacing = [0.4, 0.4, 2] +dose.hit_type = "random" + + +## SECOND DETECTOR ## +# box +# translation and rotation like in the Gate macro +box2 = sim.add_volume("Box", "box2") +box2.size = [10 * cm, 10 * cm, 105 * cm] +box2.translation = [30 * cm, 0 * cm, 52.5 * cm] +box2.material = "Vacuum" +box2.color = [0, 0, 1, 1] + +# phantoms +m = Rotation.identity().as_matrix() + +phantom2 = sim.add_volume("Box", "phantom_a_2") +phantom2.mother = "box2" +phantom2.size = [100 * mm, 100 * mm, 50 * mm] +phantom2.translation = [0 * mm, 0 * mm, -500 * mm] +phantom2.rotation = m +phantom2.material = "G4_AIR" +phantom2.color = [1, 0, 1, 1] + + +# add dose actor +dose2 = sim.add_actor("DoseActor", "doseInYZ_2") +filename = "phantom_a_2.mhd" +dose2.output = output_path / filename +dose2.mother = "phantom_a_2" +dose2.size = [250, 250, 1] +dose2.spacing = [0.4, 0.4, 2] +dose2.hit_type = "random" + +## TPS SOURCE ## +# beamline model +beamline = gate.BeamlineModel() +beamline.name = None +beamline.radiation_types = "proton" + +# polinomial coefficients +beamline.energy_mean_coeffs = [1, 0] +beamline.energy_spread_coeffs = [0.4417036946562556] +beamline.sigma_x_coeffs = [2.3335754] +beamline.theta_x_coeffs = [2.3335754e-3] +beamline.epsilon_x_coeffs = [0.00078728e-3] +beamline.sigma_y_coeffs = [1.96433431] +beamline.theta_y_coeffs = [0.00079118e-3] +beamline.epsilon_y_coeffs = [0.00249161e-3] + +# tps +nSim = 60000 # particles to simulate per beam +spots, ntot, energies, G = gate.spots_info_from_txt( + ref_path / "TreatmentPlan2Spots.txt", "proton" +) +tps = gate.TreatmentPlanSource("test", sim) +tps.set_beamline_model(beamline) +tps.set_particles_to_simulate(nSim) +tps.set_spots(spots) +tps.rotation = Rotation.from_euler("x", 90, degrees=True) +tps.initialize_tpsource() + +# add stat actor +s = sim.add_actor("SimulationStatisticsActor", "Stats") +s.track_types_flag = True + +# physics +p = sim.get_physics_user_info() +p.physics_list_name = "FTFP_INCLXX_EMZ" +sim.set_cut("world", "all", 1000 * km) +# sim.set_user_limits("phantom_a_2","max_step_size",1,['proton']) + +# create output dir, if it doesn't exist +if not os.path.isdir(output_path): + os.mkdir(output_path) + +# start simulation +sim.run() +output = sim.output + +# print results at the end +stat = output.get_actor("Stats") +print(stat) + + +# ---------------------------------------------------------------------------------------------------------------- +# tests + +# energy deposition: we expect the edep from source two +# to be double the one of source one + +print("Compare tps Edep to single pb sources") +print(" --------------------------------------- ") +mhd_1 = "phantom_a_1.mhd" +mhd_2 = "phantom_a_2.mhd" +test = True + +# check first spot +test = ( + gate.assert_images( + output_path / mhd_1, + ref_path / mhd_1, + stat, + tolerance=70, + ignore_value=0, + ) + and test +) + +# check second spot +test = ( + gate.assert_images( + output_path / mhd_1, + ref_path / mhd_1, + stat, + tolerance=70, + ignore_value=0, + ) + and test +) +print(" --------------------------------------- ") +# fig1 = gate.create_2D_Edep_colorMap(output_path / mhd_1, show=True) +# fig2 = gate.create_2D_Edep_colorMap(output_path / mhd_2, show=True) + +print("Compare ratio of the two spots with expected ratio") + +# Total Edep +is_ok = ( + gate.test_weights( + 2, + output_path / mhd_1, + output_path / mhd_2, + thresh=0.2, + ) + and test +) + + +gate.test_ok(is_ok)