diff --git a/opengate/contrib/carm/anode_heel_effect_model.py b/opengate/contrib/carm/anode_heel_effect_model.py new file mode 100644 index 000000000..a07fe5121 --- /dev/null +++ b/opengate/contrib/carm/anode_heel_effect_model.py @@ -0,0 +1,49 @@ +import numpy as np +import pandas as pd +import pathlib +import joblib +from scipy.interpolate import interp1d + + +# Function to generate phi_weights and phi_angles based on kvp (energy) +def get_phi_distribution(kvp): + # Load the scaler and model + current_path = pathlib.Path(__file__).parent.resolve() + scaler = joblib.load(current_path / "gbr_scaler.pkl") + gbr_model = joblib.load(current_path / "gbr_model.pkl") + + distances = np.arange(-9, 10, 1) + phi_angles = np.arctan(distances / 70) + + # Create a DataFrame for prediction + input_data = pd.DataFrame( + { + "Distance": distances, + "Energy": kvp, + "Theta": phi_angles, + "ThetaDegrees": np.degrees(phi_angles), + } + ) + + # Scale the input data + input_data_scaled = scaler.transform(input_data) + + # Predict the weights + phi_weights = gbr_model.predict(input_data_scaled) + + # Now interpolate between the predicted values with finer steps + fine_distances = np.arange(-9, 9.01, 0.1) + fine_phi_angles = np.arctan(fine_distances / 70) + + # Interpolate the predicted weights using linear interpolation + interp_function = interp1d(distances, phi_weights, kind="linear") + + # Interpolate predictions for the finer distances + fine_phi_weights = interp_function(fine_distances) + fine_phi_weights = fine_phi_weights[1:] + + # Adjust angles (shift them by 90 degrees in radians) + for i in range(len(fine_phi_angles)): + fine_phi_angles[i] = fine_phi_angles[i] + 180 * np.pi / 180 + + return fine_phi_weights, fine_phi_angles diff --git a/opengate/contrib/carm/anodeheeleffect.npz b/opengate/contrib/carm/anodeheeleffect.npz deleted file mode 100644 index cc60bb367..000000000 Binary files a/opengate/contrib/carm/anodeheeleffect.npz and /dev/null differ diff --git a/opengate/contrib/carm/experimental_values.json b/opengate/contrib/carm/experimental_values.json new file mode 100644 index 000000000..95baacf50 --- /dev/null +++ b/opengate/contrib/carm/experimental_values.json @@ -0,0 +1,386 @@ +{ + "50": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.0374836173001311, + 0.64720600500417, + 0.721077088049565, + 0.782199451924223, + 0.839985702371023, + 0.866198022161325, + 0.903371857500298, + 0.942213749553199, + 0.890980579053974, + 0.994042654593113, + 1.0, + 1.02823781722864, + 1.0296675801263, + 1.02311450017872, + 1.04646729417372, + 1.03991421422614, + 1.03812701060408, + 1.04253544620517, + 1.08804956511379, + 1.07887525318718, + 0.0378887167877994 + ] + }, + "60": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.0374768661320173, + 0.670573719925972, + 0.751079580505861, + 0.795496607032696, + 0.862122146822949, + 0.906539173349784, + 0.912091301665638, + 0.942319555829735, + 0.969463294262801, + 0.959901295496607, + 1.0, + 0.970080197409007, + 1.0308451573103, + 1.0410240592227, + 1.0407156076496, + 1.04380012338063, + 1.03917334978408, + 1.0613818630475, + 1.01357186921653, + 0.898519432449105, + 0.0378470080197409 + ] + }, + "70": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.0465753424657534, + 0.629144331943617, + 0.756204089735954, + 0.79630732578916, + 0.840182648401826, + 0.881080007941235, + 0.912249354774668, + 0.937661306333135, + 0.965455628350208, + 0.984117530275958, + 1.0, + 1.01290450665078, + 1.03434584077824, + 1.03911058169545, + 1.04625769307127, + 1.04744887830058, + 1.03911058169545, + 1.01111772880683, + 1.01608100059559, + 0.808616239825293, + 0.0480643240023824 + ] + }, + "80": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.0441979872663791, + 0.647360854384884, + 0.75374820291641, + 0.786609159991785, + 0.83651673855001, + 0.873279934278086, + 0.90798932018895, + 0.945163277880468, + 0.957896898747176, + 0.951735469295543, + 1.0, + 1.0145820497022, + 1.03388786198398, + 1.04066543438078, + 1.03470938591087, + 1.04805914972274, + 1.04210310125282, + 1.02752105155063, + 1.01109057301294, + 1.0262887656603, + 0.046991168617786 + ] + }, + "90": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.060656574247321, + 0.635652321823439, + 0.747576118387481, + 0.795713556727335, + 0.84504167375404, + 0.887906106480694, + 0.920734818846743, + 0.946419459091682, + 0.967171287633951, + 0.969892838918183, + 1.0, + 1.00408232692635, + 1.04643646878721, + 1.05290015308726, + 1.06072461302943, + 1.0580030617452, + 1.05409083177411, + 1.05443102568464, + 0.986052049668311, + 1.09627487667971, + 0.0656574247320973 + ] + }, + "100": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.0898376184032476, + 0.644925575101489, + 0.719891745602165, + 0.790527740189445, + 0.843031123139378, + 0.887821380243572, + 0.922733423545332, + 0.951556156968877, + 0.979296346414073, + 0.990933694181326, + 1.0, + 1.04208389715832, + 1.0553450608931, + 1.06887686062246, + 1.07510148849797, + 1.07861975642761, + 1.07307171853857, + 1.06657645466847, + 1.05439783491204, + 0.985926928281461, + 0.102692828146143 + ] + }, + "110": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.117010596765198, + 0.596319018404908, + 0.680981595092024, + 0.765755716675962, + 0.81840490797546, + 0.857557166759621, + 0.896151701059676, + 0.925153374233129, + 0.955493586168433, + 0.971556051310653, + 1.0, + 1.0228667038483, + 1.02163970998327, + 1.04138315672058, + 1.05142219743447, + 1.05543781372002, + 1.0532069157836, + 1.0499721137758, + 1.01940881204685, + 0.815393195761294, + 0.144562186279978 + ] + }, + "120": { + "distance": [ + -10, + -9, + -8, + -7, + -6, + -5, + -4, + -3, + -2, + -1, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10 + ], + "weight": [ + 0.119681149657573, + 0.617716402829235, + 0.691366341080049, + 0.77500842034355, + 0.82519366790165, + 0.867070843157067, + 0.904120354777141, + 0.934770405299203, + 0.960368249691254, + 0.975861681823285, + 1.0, + 1.02234197821938, + 1.04064219153475, + 1.05029751880543, + 1.05883013360278, + 1.06264735601213, + 1.05838104861345, + 1.0511956887841, + 1.02627147187605, + 0.771977096665544, + 0.145728079038958 + ] + } +} \ No newline at end of file diff --git a/opengate/contrib/carm/gbr_model.pkl b/opengate/contrib/carm/gbr_model.pkl new file mode 100644 index 000000000..d83f0a25e Binary files /dev/null and b/opengate/contrib/carm/gbr_model.pkl differ diff --git a/opengate/contrib/carm/gbr_scaler.pkl b/opengate/contrib/carm/gbr_scaler.pkl new file mode 100644 index 000000000..3f58e4edb Binary files /dev/null and b/opengate/contrib/carm/gbr_scaler.pkl differ diff --git a/opengate/contrib/carm/siemensciosalpha.py b/opengate/contrib/carm/siemensciosalpha.py index f45e83bf7..01f987af3 100644 --- a/opengate/contrib/carm/siemensciosalpha.py +++ b/opengate/contrib/carm/siemensciosalpha.py @@ -2,6 +2,7 @@ from scipy.spatial.transform import Rotation import numpy as np from opengate.utility import LazyModuleLoader +from .anode_heel_effect_model import get_phi_distribution sp = LazyModuleLoader("spekpy") import pathlib @@ -111,12 +112,11 @@ def add_carm_source(self, kvp): source.direction_relative_to_attached_volume = True source.direction.type = "histogram" source.direction.histogram_theta_weights = [1] - source.direction.histogram_theta_angles = [85 * deg, 95 * deg] + source.direction.histogram_theta_angles = [82.68 * deg, 97.32 * deg] - # TODO: Need real values for the anode heel effect - data = np.load(current_path / "anodeheeleffect.npz") - source.direction.histogram_phi_weights = data["weight"] - source.direction.histogram_phi_angles = data["angle"] + phi_weights, phi_angles = get_phi_distribution(kvp) + source.direction.histogram_phi_weights = phi_weights + source.direction.histogram_phi_angles = phi_angles source.energy.type = "histogram" source.energy.histogram_weight = weights @@ -137,20 +137,20 @@ def add_collimators(self): collimators = [ { - "translation": [37.5 * mm, 0 * cm, -z_xray_tank / 2 * mm + 1 * mm], - "size": [2.5 * cm, 5 * cm, 1 * mm], + "translation": [42.5 * mm, 0 * cm, -z_xray_tank / 2 * mm + 1 * mm], + "size": [1.5 * cm, 10 * cm, 1 * mm], }, { "translation": [-37.5 * mm, 0 * cm, -z_xray_tank / 2 * mm + 1 * mm], - "size": [2.5 * cm, 5 * cm, 1 * mm], + "size": [1.5 * cm, 10 * cm, 1 * mm], }, { "translation": [0 * cm, 37.5 * mm, -z_xray_tank / 2 * mm + 3 * mm], - "size": [5 * cm, 2.5 * cm, 1 * mm], + "size": [10 * cm, 1.5 * cm, 1 * mm], }, { "translation": [0 * cm, -37.5 * mm, -z_xray_tank / 2 * mm + 3 * mm], - "size": [5 * cm, 2.5 * cm, 1 * mm], + "size": [10 * cm, 1.5 * cm, 1 * mm], }, ] @@ -167,7 +167,7 @@ def add_collimators(self): killer.attached_to = [f"{self.machine_name}_collimator{i+1}" for i in range(4)] def set_collimation(self, collimation1, collimation2): - if not 0 <= collimation1 <= 25 or not 0 <= collimation2 <= 25: + if not 10 <= collimation1 <= 25 or not 10 <= collimation2 <= 25: raise ValueError("Collimation values must be between 0 and 25 mm") collimation1 = 25 - collimation1 @@ -177,10 +177,10 @@ def set_collimation(self, collimation1, collimation2): z_xray_tank = xray_tank.size[2] translations = [ - [37.5 * mm - collimation1, 0 * cm, -z_xray_tank / 2 * mm + 1 * mm], - [-37.5 * mm + collimation1, 0 * cm, -z_xray_tank / 2 * mm + 1 * mm], - [0 * cm, 37.5 * mm - collimation2, -z_xray_tank / 2 * mm + 3 * mm], - [0 * cm, -37.5 * mm + collimation2, -z_xray_tank / 2 * mm + 3 * mm], + [42.5 * mm - collimation1, 0 * cm, -z_xray_tank / 2 * mm + 1 * mm], + [-42.5 * mm + collimation1, 0 * cm, -z_xray_tank / 2 * mm + 1 * mm], + [0 * cm, 42.5 * mm - collimation2, -z_xray_tank / 2 * mm + 3 * mm], + [0 * cm, -42.5 * mm + collimation2, -z_xray_tank / 2 * mm + 3 * mm], ] for i, translation in enumerate(translations): diff --git a/opengate/tests/src/geometry/test075_siemens_cios_alpha.py b/opengate/tests/src/geometry/test075_siemens_cios_alpha.py index 666c9d6d6..93caadb1d 100755 --- a/opengate/tests/src/geometry/test075_siemens_cios_alpha.py +++ b/opengate/tests/src/geometry/test075_siemens_cios_alpha.py @@ -4,22 +4,31 @@ import opengate as gate from opengate.contrib.carm.siemensciosalpha import Ciosalpha from scipy.spatial.transform import Rotation +from opengate.tests import utility as utl +import matplotlib.pyplot as plt +import itk +import numpy as np +import json +import pathlib if __name__ == "__main__": # paths - # paths = utility.get_default_test_paths(__file__, output_folder="test075_siemens_cios_alpha") + paths = utl.get_default_test_paths( + __file__, "gate_test075_siemens_cios_alpha", "test075" + ) # create the simulation sim = gate.Simulation() # main options sim.g4_verbose = False - # sim.visu = True + sim.visu = False sim.visu_type = "vrml" sim.number_of_threads = 1 sim.random_seed = 12345678 sim.check_volumes_overlap = True + sim.output_dir = paths.gate_output # units m = gate.g4_units.m @@ -33,27 +42,153 @@ world.material = "G4_AIR" # xray tube spectrum parameters - # tube potential [kV] - kvp = 100 + # tube potential [kV] (50, 60, 70, 80, 90, 100, 110, 120 to compare with experimental data) + kvp = 80 # add a carm carm = Ciosalpha(sim, kvp, source_only=False) - carm.rotation = Rotation.from_euler("ZYX", [0, 20, 0], degrees=True).as_matrix() + carm.rotation = Rotation.from_euler("ZYX", [0, 90, 0], degrees=True).as_matrix() carm.translation = [0 * cm, 0 * cm, 0 * cm] carm.collimation = [25 * mm, 25 * mm] - carm.source.n = 1e6 + carm.source.n = 1e7 if sim.visu: carm.source.n = 1000 - # aluminum table - table = sim.add_volume("Box", "aluminum_table") - table.size = [60 * cm, 2 * m, 0.9 * mm] - table.material = "G4_Al" - table.translation = [0 * m, 0 * cm, 0 * cm] - table.color = [0.8, 0.8, 0.8, 1] + # CBCT detector plane + detector_plane = sim.add_volume("Box", "CBCT_detector_plane") + detector_plane.size = [10 * mm, 200 * mm, 200 * mm] + detector_plane.translation = [-150 * mm, 0, 0] + detector_plane.material = "G4_AIR" + detector_plane.color = [1, 0, 0, 1] + + # physics + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option1" + sim.physics_manager.set_production_cut("world", "all", 10 * mm) + + # actor + detector_actor = sim.add_actor("FluenceActor", "detector_actor") + detector_actor.attached_to = detector_plane + detector_actor.output_filename = "fluence.mhd" + detector_actor.spacing = [10 * mm, 2 * mm, 2 * mm] + detector_actor.size = [10, 101, 101] + detector_actor.output_coordinate_system = "local" # start simulation sim.run() - # TODO: Test + def get_normalized_profile(img, axis="z", nb_points=200): + """ + Get evenly spaced values from the profile and normalize them by the middle value. + + Parameters: + ----------- + img : itk.Image + The input image + axis : str + The axis along which to get the profile ('x', 'y', or 'z') + nb_points : int + Number of points to sample (default=200) + + Returns: + -------- + tuple: (positions, normalized_values) + The positions and normalized values + """ + # Get data in numpy array + data = itk.GetArrayViewFromImage(img) + + if axis == "z": + y = np.nansum(data, 2) + y = np.nansum(y, 1) + x = np.arange(len(y)) * img.GetSpacing()[2] + elif axis == "y": + y = np.nansum(data, 2) + y = np.nansum(y, 0) + x = np.arange(len(y)) * img.GetSpacing()[1] + else: # axis == "x" + y = np.nansum(data, 1) + y = np.nansum(y, 0) + x = np.arange(len(y)) * img.GetSpacing()[0] + + # Get middle value for normalization + mid_idx = len(x) // 2 + mid_val = y[mid_idx] + + # Interpolate to get evenly spaced points + from scipy.interpolate import interp1d + + f = interp1d(x, y, kind="linear") + + # Create new x points evenly spaced + x_new = np.linspace(x[0], x[-1], nb_points) + y_new = f(x_new) + + # Normalize by middle value + y_normalized = y_new / mid_val + + return x_new, y_normalized + + # Load experimental and simulation data + # Simulation data + img_mhd_out = itk.imread(paths.gate_output / detector_actor.output_filename) + positions, weights = get_normalized_profile(img_mhd_out, axis="z", nb_points=200) + positions = [pos - 100 for pos in positions] + + # Experimental data + current_dir = pathlib.Path(__file__).parent + exp_file_path = ( + current_dir.parent.parent / "contrib" / "carm" / "experimental_values.json" + ) + with open(exp_file_path, "r") as f: + exp_data = json.load(f) + exp_positions = [pos * 10.0 for pos in exp_data[str(kvp)]["distance"]] + exp_weights = list(reversed(exp_data[str(kvp)]["weight"])) + + # Calculate differences + differences = [] + within_uncertainty = [] + print(f"\nAnalysis for {kvp} kV:") + print("Position (mm) | Simulation | Experimental | Difference (%)") + print("-" * 75) + + # Interpolate simulation data to match experimental positions + from scipy.interpolate import interp1d + + f = interp1d(positions, weights, kind="linear") + sim_weights_interp = f(exp_positions) + + # excluding first/last 2 values + exp_positions_trimmed = exp_positions[2:-2] + sim_weights_interp_trimmed = sim_weights_interp[2:-2] + exp_weights_trimmed = exp_weights[2:-2] + + for pos, sim_w, exp_w in zip( + exp_positions_trimmed, sim_weights_interp_trimmed, exp_weights_trimmed + ): + diff_percent = ((sim_w - exp_w) / exp_w) * 100 + differences.append(diff_percent) + print(f"{pos:13.1f} | {sim_w:10.4f} | {exp_w:12.4f} | {diff_percent:11.2f}%") + + # results + max_difference = max(abs(d) for d in differences) + avg_abs_difference = sum(abs(d) for d in differences) / len(differences) + print(f"\nMaximum difference: {max_difference:.2f}%") + print(f"Average absolute difference: {avg_abs_difference:.2f}%") + + # Create plot with updated title + fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(25, 10)) + ax.plot(positions, weights, "b-", label="Simulation") + ax.plot(exp_positions, exp_weights, "ro", label=f"Experimental ({kvp} kV)") + ax.set_xlabel("Position (mm)") + ax.set_ylabel("Weight (normalized to middle value)") + ax.set_title(f"Anode Heel Effect Profile Comparison - Siemens Cios Alpha {kvp} kV") + ax.grid(True) + ax.legend() + + # Save plot + fig.savefig(paths.gate_output / "anode_heel_effect_comparison.png") + plt.show() + + is_ok = max_difference < 5.0 and avg_abs_difference < 2.0 + utl.test_ok(is_ok)